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Abstract 


GPS PPP has found many scientific and industrial applications due to its cost- 
effectiveness, global coverage, and high accuracy for decades. However, it suffers from 
a few drawbacks which limites further applications, e.g., slow convergence and 
unresolved ambiguity. The rapid development of multi-GNSS and multi-frequency 
signals offers exciting prospects for further improvements. In this thesis, multi- 
frequency and multi-GNSS measurements are properly integrated for PPP with 
ambiguity resolution to achieve the state-of-the-art fast and accurate positioning, which 
provides an important contribution to GNSS precise positioning and applications. 


Integer ambiguity resolution is the key issue for improving PPP accuracy and 
convergence. The core of PPP ambiguity resolution is carrier phase fractional cycle biases 
(FCB) estimation. In this thesis, the characteristics of phase FCB are analyzed with the 
recent and long temporal series. Results indicate that the temporal stabilities have been 
significantly improved in recent years. Taking advantage of this property, an improved 
FCB estimation method based on Kalman filter is proposed. The designed Kalman filter 
significantly reduces the dimension of the involved matrix and accelerates FCB 
computation, which outperforms the commonly used least square method in terms of 
efficiency. The method is especially useful for real-time application. 


With the estimated FCB products, the PPP ambiguity resolution with the current Galileo 
and BDS constellation is verified. The satellite FCB is thoroughly assessed by a 
comparison with that of GPS in terms of data usage rate, residual distribution, as well as 
standard deviation of daily estimates. Results indicate that the quality of Galileo wide- 
lane (WL) FCB is better than GPS and BDS, while the quality of Galileo narrow-lane (NL) 
FCB is slightly worse than that of GPS but better than that of BDS. Within the Galileo 
constellation, the performance of In-Orbit Validation (IOV) satellites WL FCB is worse 
than that of Full Operational Capability (FOC) satellites as a result of a reduction in the 
power of the transmitted signal. The WL FCB performance of the two highly eccentric 
satellites is comparable to other FOC satellites. The NL FCB quality of FOC, IOV (except 
E19), as well as the two eccentric satellites, shows no significant difference in terms of 
data usage rates and residuals. Solution of PPP ambiguity resolution demonstrates that 
the Galileo and BDS observations can bring an obvious benefit to GPS only PPP 
ambiguity resolution. 


PPP ambiguity resolution is also extended to multi-frequency observations in this thesis. 
In order to properly integrate the multi-frequency observations, a unified uncombined 
PPP ambiguity resolution model based on raw observations is proposed. Based on the 
unified model, the FCBs generated from multi-frequency observations can be flexibly 
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used, such as for dual- or triple- frequency PPP AR. Its efficiency is verified with Galileo 
and BeiDou triple-frequency observations collected from globally distributed MGEX 
stations. The estimated FCB are assessed with respect to residual distributions and 
standard deviations, which shows a good consistency with the input float ambiguities. 
The standard deviation of linear combined FCB is much smaller than that of raw FCB, 
which indicates linear combined FCB is more efficient for broadcasting in real time 
applications. Triple-frequency PPP ambiguity resolution indicates that the positional 
biases are significantly reduced compared with the float solutions. The improvements 
are 49.2%, 38.3%, and 29.6%, respectively, in east/north/up components for positioning 
with BDS, while the corresponding improvements are 60.0%, 29.0%, and 21.1% for 
positioning with Galileo. 


The uncombined signal processing brings new challenges for cycle slip detection. The 
most significant feature is carrier frequency identification of cycle slips. Since all carrier 
frequency observations are processed separately, it is essential to identify the carrier 
frequency of the cycle slips to avoid contaminating other observations. To provide 
continuous carrier phase measurements for the above research, an improved approach 
based on a time-differenced model for cycle slip detection and repair is proposed, which 
reduce the false alarms and increase the success rate of cycle slip estimation. 


In summary, the GPS only PPP has been extended to multi-frequency and multi-GNSS 
PPP ambiguity resolution with improved accuracy and fast convergence. This is 
accomplished by a unified model based on the uncombined PPP. The proposed model 
has been carefully studied and enriched with improved cycle slip detection and repair, 
fast FCB estimation method and ambiguity resolution. 


This article-based (cumulative) thesis consists of a detailed introductory chapter and 
four chapters associated with the following peer-reviewed publications: 


1) Xiao G, Li P, Gao Y, Heck B (2019a) A Unified Model for Multi-Frequency PPP 
Ambiguity Resolution and Test Results with Galileo and BeiDou Triple- 
Frequency Observations. Remote Sensing 11(2):116 doi:10.3390/rs11020116 


2) Xiao G, Li P, Sui L, Heck B, Schuh H (2019b) Estimating and assessing Galileo 
satellite fractional cycle bias for PPP ambiguity resolution. GPS Solutions 23:3 
doi:10.1007/s10291-018-0793-z 


3) Xiao G, Sui L, Heck B, Zeng T, Tian Y (2018b) Estimating satellite phase fractional 
cycle biases based on Kalman filter. GPS Solutions 22:82 doi:10.1007/s10291-018- 
0749-3 


4) Xiao G, Mayer M, Heck B, Sui L, Zeng T, Zhao D (2018a) Improved time- 
differenced cycle slip detect and repair for GNSS undifferenced observations. 
GPS Solutions 22:6 doi:10.1007/s10291-017-0677-7 


Zusammenfassung 


Aufgrund seiner Wirtschaftlichkeit, globalen Abdeckung und hohen Genauigkeit hat 
GPS PPP seit Jahrzehnten. PPP weist jedoch einige Nachteile auf, die weitere 
Anwendungen einschränken, z. B. langsame Konvergenzzeit und ungelöste Mehr- 
deutigkeit. Die schnelle Entwicklung von GNSS und seinen Mehrfrequenz-Signalen 
bietet spannende Perspektiven fiir weitere Verbesserungen. In dieser Dissertation sind 
Multi-GNSS- und Mehrfrequenz-Beobachtungen für PPP-Lösung mit Mehrdeutig- 
keitauflosung integriert, um eine schnelle und genauige Positionsbestimmung auf dem 
neuesten Stand der Technik zu erreichen. Dies ist ein wichtiger Beitrag zur präzisen 
Positionsbestimmung mittels GNSS und seinen Anwendungen. 


Festsetzung der Phasenmehrdeutigkeit auf ihre ganzzahlige Werte ist von großer 
Bedeutung fiir die Verbesserung der PPP-Genauigkeit und -Kovergenzzeit. Der 
Schwerpunkt in diesem Zusammenhang ist die Schätzung der (Fractional Cycle Biases) 
FCB der Tragerphasen. In dieser Dissertation-werden die Eigenschaften der Phasen-FCB 
mit kürzlichen und langzeitlichen Serien analysiert. Die Ergebnisse zeigen, dass die 
zeitliche Stabilitat in den letzten Jahren wesentlich verbessert wird. Bei Nutzung dieser 
Eigenschaft wird ein verbessertes FCB-Schatzverfahren mit Kalman-Filter entwickelt. 
Dieser entwicklende Kalman-Filter reduziert die Dimension der betroffenen Matrix und 
beschleunigt damit die Berechung des FCB. Allerdings verliert die traditionale LS 
Methode hinsichtlich der Effizienz im Vergleich dazu. diese Methdode eignet sich 
besondersfür die Echtzeitiganwendungen. 


Mit den geschätzten FCB-Produkten wird die PPP-Mehrdeutigkeitauflösung mit der 
aktuellen Konstellation von Galileo und BDS überprüft. DerStatelliten-FCB wird im 
Vergleich mit GPS gewonnen Wertenhinsichtlich der Datennutzungsrate, der Verteilung 
der Residuen und der taglichen Abschatzungen der Standardabweichung bewertet. Den 
Ergebnisse zeigen, dass die Qualitat der Galileo Wide-Lane (WL) FCB besser als GPS 
und BDS ist. Mittlerweile ist die Qualitat der Galileo Narrow-Lane (NL) FCB ein bisschen 
schlechter als GPS aber besser als BDS. innerhalb der Galileo Konstellation sind die 
Leistung des In-Orbit-Validation (IOV)-Satelliten bezüglich WL FCB schlechter als die 
Full Operational Capability (FOC) Stalliten wegen der Dampfung der Energie von den 
austrahlenden Signalen. Die WL FCB Leistung der Satelliten mit hohem Exzentrizitat ist 
vergleichbar mit den anderen FOC-Satelliten. Wenn die Datennutzungsrate und 
Residuen betrachtet werden, weist die NL FCB Qualität von FOC, IOV (Außer E19), und 
den anderen Statelliten mit hohem Exzentrizitat keine wesentliche Differenz auf. Die 
PPP Mehrdeutigkeitauflösung zeigt, dass die Galileo und BDS Beobachtungen einen 
deutlichen Vorteil für GPS-only PPP Mehrdeutigkeitauflösung darstellen können. 


In dieser Dissertation wird die PPP-Mehrdeutigkeitsauflösung auch auf Mehr- 
frequenzbeobachtungen erweitert. Um die Mehrfrequenzbeobachtungen richtig zu 
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integrieren, wird basierend auf Rohbeobachtungen ein einheitliches unkombiniertes 
PPP-Mehrdeutigkeitsauflösungsmodell entwickelt. Mit solchem einheitlichen Modell 
können die von Mehrfrequenz-Beobachtungen erzeugten FCBs flexibel verwendet 
werden, beispielsweise für zwei- oder drei Frequenzen PPP Mehrdeutigkeitauflösung. 
Die Effizienz dieses Modells wird mit Galileo und Beidou drei-Frequenz-Beobachtungen 
geprüft, die von globalen verteilten MGEX Stationen gesammelt werden. Die 
geschätzten FCBs werden hinsichtlich der Residuenausverteilung und Standard- 
abweichungen bewertet, was eine gute Übereinstimmung mit der float-Lösung der 
Mehrdeutigkeit zeigt. Die Standardabweichung von linear kombinierten FCBs ist viel 
geringer als die von Roh-FCBs. Diese bedeutet, dass ein linear kombinierte FCBs 
effizienter für das Broadcasting in Echtzeitanwendungen sind. Die drei-Frequenz-PPP- 
Mehrdeutigkeitauflösung zeigt, dass die Positionsabweichung im Vergleich zur float- 
Lösung stark reduziert sind. Die Verbesserungen für die Ost, Nord und Höhen- 
kompenenten mit BDS Positionierung sind 49.2%, 38.3% und 29.6%, in der Reheinfolge. 
Die entsprechenden Verbesserungen mit Galileo Positionierung sind 60.0%, 29.0% und 
21.1%. 


Neue Herausforderung für Cycle-slip Erkennung entsteht durch die Verarbeitung des 
unkombinierten Signals. Da alle Tragerfrequenzbeobachtungen separat verarbeitet 
werden, ist es sehr wichtig, die Trägerfrequenz des Cycle slips zu identifizieren, damit 
die anderen Beobachtungen nicht davon beeinflusst werden können. Um konti- 
nuierliche Trägerphasenmessungen für die vorher genannten Forschungen bereit- 
zustellen, wird ein verbesserter Ansatz entwickelt, der auf einem zeitlich differenzierten 
Modell für die Erkennung und Reparatur des Cycle-slips basiert. Dieses Modell 
reduziert die Gefahr für falsche Entscheidung und erhöht die Chance auf eine 
erfolgreiche Abschätzung des Cycle-slips. 


Zusammenfassend ist festzuhalten, dass die PPP mittles only-GPS auf Mehrfrequenz 
und Mehr-GNSS-PPP-Mehrdeutigkeitauflösung mit verbesserter Genauigkeit und 
schneller Konvergenz erweitert wurde. Dies wurde durch ein einheitlichtes Modell 
erreicht, das auf nicht kombinierte PPP basiert. Dieses Model wurde sorgfältig 
untersucht und mit verbesserter Erkennung und Reparatur von cycle-slip, schneller 
FCB-Schätzverfahren und Mehrdeutigkeitsauflösung angereichert. 


Diese artikel-basierte (kumulative) Dissertation besteht aus einem ausführlichen 
Einführungskapitel und vier Kapiteln, die mit den folgenden Publikationen verbunden 
sind: 
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1) Xiao G, Li P, Gao Y, Heck B (2019a) A Unified Model for Multi-Frequency PPP 
Ambiguity Resolution and Test Results with Galileo and BeiDou Triple- 
Frequency Observations. Remote Sensing 11(2):116 doi:10.3390/rs11020116 


2) Xiao G, Li P, Sui L, Heck B, Schuh H (2019b) Estimating and assessing Galileo 
satellite fractional cycle bias for PPP ambiguity resolution. GPS Solutions 23:3 
doi:10.1007/s10291-018-0793-z 


3) Xiao G, Sui L, Heck B, Zeng T, Tian Y (2018b) Estimating satellite phase fractional 
cycle biases based on Kalman filter. GPS Solutions 22:82 doi:10.1007/s10291-018- 
0749-3 


4) Xiao G, Mayer M, Heck B, Sui L, Zeng T, Zhao D (2018a) Improved time- 
differenced cycle slip detect and repair for GNSS undifferenced observations. 
GPS Solutions 22:6 doi:10.1007/s10291-017-0677-7 
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Chapter 1 
Introductory chapter 


1 Introduction 


1.1 Motivation and background 


Satellite-based positioning has become an integral part of our modern-day society and 
is used by people all over the world (Teunissen and Montenbruck 2017). For positioning 
with satellite systems, pseudorange and carrier phase are the basic observables. Due to 
the low precision of pseudorange measurements and broadcast ephemeris, the resulting 
accuracy of standard point positioning is at several meters level (Hofmann-Wellenhof et 
al. 2007). For high-accuracy positioning at centimeters level, the carrier phase measure- 
ments must be adopted. The carrier phase-based positioning techniques can be roughly 
divided as differential positioning, such as real time kinematic (RTK), and absolute po- 
sitioning, i.e., precise point positioning (PPP). Differential positioning requires reference 
receivers to mitigate the common errors, while PPP is characterized by single receiver 
with sophisticated error models and precise satellite information. PPP is able to achieve 
centimeter to sub centimeters accuracies with dual-frequency carrier phase and pseu- 
dorange measurements. Due to its cost-effectiveness, global coverage, and high accu- 
racy, PPP has found increased applications. However, PPP shows a few drawbacks 
(Bisnath and Gao 2008), e.g., slow convergence and lower accuracy compared to differ- 
ential positioning. The integer ambiguity resolution (AR) technique is expected to fur- 
ther enhance the accuracy and shorten the convergence time. Although the theory was 
proposed and elaborated in (Gabor and Nerem 1999; Gao and Shen 2002), PPP AR was 
not practically realized and implemented until (Collins et al. 2008; Ge et al. 2008; 
Laurichesse et al. 2009). It is proved that AR is able to improve the PPP accuracy, espe- 
cially the east component, to a comparable level as that of differential positioning. Many 
applications including ionospheric modelling (Banville and Langley 2011), tropospheric 
modelling (Li et al. 2015c), time transfer (Petit et al. 2015), seismic displacement moni- 
toring (Li et al. 2013c), low earth satellite orbit (Montenbruck et al. 2018) and navigation 
satellite orbit determination (Katsigianni et al. 2018), have further confirmed its effi- 
ciency. After ten years development, the technique is so impressive that the global navi- 
gation satellite system (GNSS) community proposed to establish a specific working 
group for PPP AR at the IGS Workshop 2018. 


While remarkable progress has been made with legacy GPS observations, the ongoing 
modernization and the build-up of new GNSS constellations (shown in Table 1) offers 
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exciting prospects for further improvements. The GNSS has evolved from GPS to four 
major global systems including the Russian GLONASS, Chinese BDS, and European 
Galileo. Specifically, GLONASS has restored its full orbital constellation of 24 satellites 
since 2011 (Revnivykh 2012). For BDS, the construction is planned for three steps (Yang 
et al. 2018). The formal establishment of BDS-1 was marked by the launch of the third 
BeiDou satellite in 2003. The second step is the BDS-2 regional system, which was ac- 
complished by a constellation of 14 satellites in 2012. With the successful launch of BDS 
third generation experimental satellites in 2015, the construction of the BDS global sys- 
tem (BDS-3) starts to speed up, and is expected to be completed in 2020. Similarly, Gali- 
leo also witnesses remarkable achievements in recent years (Steigenberger and 
Montenbruck 2017). On 15 December 2016, Galileo started offering initial services. As of 
March 2019, there are 22 satellites in orbit and 4 satellites under commissioning. The 
complete constellation with 30 satellite could be expected by 2020. 


As can be seen from Table 1, there are 110 GNSS satellites in orbit with 10 satellites under 
commissioning at present. When all four global systems reach their full constellations, 
there will be more than 120 usable satellites. Compared with 32 GPS satellites, the fusion 
of multiple GNSS can significantly increase the number of visible satellites by a factor of 
nearly 4 and improve the spatial geometry. Fast convergence, improved accuracy and 
enhanced reliability can be expected. Under this circumstance, the contribution of multi- 
GNSS integration to GPS standalone PPP, as well as ambiguity resolution, should be 
investigated. 


Table 1 Status of multi-GNSS as of Mar. 2019 


System GPS GLONASS Galileo BDS 
No of üsable 31 94 32 33 
satellites 
No of satellite 
in commission- N/A N/A 4 6 
ing 
Full constella- 
tion and ex- 24, Already 24, Already 30, 2020 35,2020 
pected date 
Early Operational Capa- 


On 27 Decem- 
ber 2018, BDS 


Launch of the The K series, which 
first GPS III support CDMA sig- 
Remark satellites on 23 nals, should gradu- 


bility (EOC) on 15 De- 


cember 2016, and is ex- 
started to pro- 


ted t h Full 
December ally replace existing pec A ae b vide global ser- 
2018 satellite en pa vices 
(FOC) in 2019 


In recent years, PPP has found new opportunity in multi-frequency observations. The 
additional signals can be used to form combinations with excellent properties, such as 
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low noise and reduced ionospheric delay (Feng 2008). The additional signals are ex- 
pected to further enhance the performance of PPP (Geng and Bock 2013). As of March 
2019, triple- or multi-frequency signals have become available from more and more 
GNSS satellites, as shown in Table 2. In the effort to modernize GPS, new L5 signals were 
added, and are available now on all satellites from GPS IIF (Tegedor and Ovstedal 2014). 
The L5 signals provide secure and robust transmissions for life critical applications. 
However, there exists additional inter-frequency clock bias (Montenbruck et al. 2012), 
which should be accounted for (Li et al. 2013a; Pan et al. 2018). The legacy GLONASS 
satellites transmit on a different frequency using a 15-channel frequency division multi- 
ple access (FDMA) technique (Glonass 2008). The 24-satellite constellation is accommo- 
dated with only 15 channels, as identical frequency channels can be used for antipodal 
satellite pairs which are never both in view at the same time for earth-based users. The 
FDMA technique introduces the receiver inter-frequency bias (IFB) between satellites, 
which is a problem for high-precision positioning (Reußner and Wanninger 2012; 
Wanninger 2012). To avoid the deficiency, new code division multiple access (CDMA) 
signals were researched for GLONASS since 2008. The GLONASS-K1 satellites launched 
in 2011 introduced L3 CDMA signals, followed by the enhanced GLONASS-K1 and 
GLONASS-K2 satellites in 2018 which introduced L1, L2 and L3 CDMA signals (IAC 
2016). These CDMA signals may replace the FDMA signals in the future. Galileo pro- 
vides an open signal in the E1 band and a wideband signal covering the E5 a&b band 
(Zaminpardaz and Teunissen 2017). The Alternating Binary Offset Carrier (AltBOC) 
modulation can either be tracked as a composite signal or as distinct signals in the E5a 
and E5b sub-bands. For BDS-1 and BDS-2, three signals, e.g., B1, B2, and B3, were broad- 
casted, while three new carrier frequencies (B1C, B2a, Bs) were added in BDS-3 (Yang et 
al. 2018). 


Table 2 Status of carrier frequency of GNSS as of Mar. 2019 [Unit: MHz] 


System GPS GLONASS Galileo BDS 
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The multi-frequency signals are anticipated to enable rapid ambiguity resolution in dif- 
ferential positioning. Results indicate that instantaneous triple-frequency ambiguity res- 
olution can be achieved with optimal combinations (Jung et al. 2000; Vollath et al. 1999). 
The underlying foundation is that combinations with reduced ionospheric delays, low 
noise and longer wavelength can be formed from multi-frequency observations (Geng 
and Bock 2013). Considering the proved efficiency in differential positioning, multi-fre- 
quency observations are expected to improve the performance of PPP. 


In conclusion, multi-frequency and multi-constellation GNSS observations show a clear 
advantage over standard dual-frequency GPS observations. For PPP, improved accuracy 
and fast convergence time can be expected. However, proper integration of measure- 
ments from multi-constellation and multi-frequency GNSS into PPP requires detail in- 
vestigations, which motivate this research. 


1.2 Outline 


The outline of the rest of this introductory chapter is presented as follows: Section 2 in- 
troduces the traditional dual-frequency PPP and ambiguity resolution models. Section 3 
highlights the limitations of the current model, followed by the research objectives of 
this thesis. While Section 4 summarizes the main contribution of this thesis, Section 5 
provides an overview of the related publications as presented in Chapters 2-5. This in- 
troductory chapter closes with Section 6, in which an overall conclusion of the thesis 
including an outlook for the future research is presented. 


2 Foundations 


For a satellite s observed by receiver r, the corresponding raw pseudorange and carrier 
phase observation equations can be expressed as (Hofmann-Wellenhof et al. 2007; Leick 
et al. 2015; Xu 2007) 


|; cpi dt, at Par + ap dle, + Dag = De d, 
(1) 


nf = pe + dt, — dt’ + dT — ag: dlp, - Aj(N2g + Bp — Bf) + £o, 


where the subscript f = (1,253, --- ) refers to a specific carrier frequency, superscript s re- 
fers to a specific satellite; p? indicates the geometric distance between the satellite and 
receiver; dt, and dt? are the clock errors of receiver and satellite; dT is the slant tropo- 
spheric delay; dl; is the slant ionospheric delay on the first carrier frequency and ay = 
Af / M is the carrier frequency-dependent factor; D,,- and Dj are the receiver and satel- 
lite specific code hardware delays; Ay and N; p are the wavelength in meter and integer 
ambiguity in cycle; B. and Bj are the receiver-dependent and satellite-dependent un- 
calibrated phase delays; ep, and £5, are the pseudorange and carrier phase measure- 
ment noise, respectively. Note that the higher-order ionospheric effects are neglected, as 
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they have limited influence on the performance of ambiguity resolution (Hadas et al. 
2017). 


In the classic GPS dual-frequency PPP, the ionospheric-free (IF) combination is routinely 
employed to eliminate the effect of the first-order ionospheric delay(Zumberge et al. 
1997). After multiplied by the IF coefficients, the corresponding pseudorange and carrier 
phase observation equations can be expressed as 


TM = p; + dt, — dt + dT + D, jp — Dir + €p,, (2) 


PF ip = D? + dt, T dts + dT + Ai (Ne ir + Br ir — B.) + Edif 


Conventionally, precise orbit and clock products from the international GNSS service 
(IGS) analysis center are used to remove satellite orbit and clock errors. The pseudorange 
ionospheric-free hardware delay bias Dj; is assimilated into the clock offset dt? in ac- 
cordance with the IGS analysis convention. Due to the fact that pseudorange measure- 
ments provide the reference to clock parameters, the actual receiver clock estimate 
would absorb the ionospheric-free combination of the receiver pseudorange hardware 
delay D, r. After applying the GNSS precise satellite clock products, equation (2) can be 
rewritten as 


PS ip = ps + dT + dt, + &,, 


= = (3) 
PF ir = p? + dT + dt, + dA t Edif 
where the receiver clock and ambiguity can be re-parameterized as 
dt, = (dt, + Dyr) (4) 
XS E S 
Near = Nr ır + brir — big (5) 
Dyjg 
b, ig = Briar — ee (6) 
IF 
Dir 
bîp = Bir — Ir (7) 


IF 


With at least four satellites simultaneously observed, the equations can be integrated and 
resolved with least squares algorithm. The estimable parameters are receiver coordi- 
nates, tropospheric delay, receiver clock offset and ambiguities. Besides this estimation 
strategy, the state-of-the-art corrections (Kouba 2009; Zumberge et al. 1997) including 
the satellite and receiver antenna phase center offset and variation (Schmid et al. 2005), 
the relativistic clock effects (Kouba 2004), the phase wind-up (Wu et al. 1991), the site 
displacement effects due to solid earth tide and ocean loading (Petit and Luzum 2010), 
and the differential code biases (Dach et al. 2015) should be considered. Fortunately, 
these biases are well researched and can be corrected with models or external products. 
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The above discussed model is referred to as standard dual-frequency PPP model. Note 
that the estimated ambiguity parameter is a combination of the integer ambiguity, the 
corresponding code hardware delays, and the uncalibrated carrier phase delays at both 
receiver and satellite ends. In this case, the integer property is lost. The fact that double- 
differenced ambiguities in global or regional networks can be resolved to integer values 
(Dong and Bock 1989) lays the foundation of integer ambiguity resolution for PPP. The 
resolved double-differenced ambiguity implies that the fractional parts of two single- 
differenced ambiguities (across satellites) must agree well with each other (Ge et al. 
2008). By estimating the fractional parts, denoted as fractional cycle bias (FCB) in the 
sequel, at the server end and applying them to single differencing PPP at the user end, 
PPP integer ambiguity resolution can be achieved. 


For PPP ambiguity resolution, the term N, ir is usually decomposed into the following 
combination 


ch 


ER Fe _ fe Newt + LL 

rIF = =. ge (8) 
where N72, is the integer wide-lane (WL) ambiguity and N, uL is the float narrow-lane 
(NL) ambiguity. Usually, the WL ambiguity is resolved by the Hatch-Melbourne-Wüb- 
bena (HMW) combination observable (Hatch 1982; Melbourne 1985; Wübbena 1985). 
With the fixed WL ambiguity, the float NL ambiguity can be derived based on (8), and 
tested whether it is also fixable. An ionospheric-free ambiguity is fixed only when both 
its WL and NL ambiguities are fixed. Now, the fixing of IF ambiguity has been trans- 
formed to fixing of WL and NL ambiguities. Similarly, the estimation of IF FCB has been 
transformed to estimation of WL and NL FCBs. The method of WL and NL FCB deter- 
minations is presented in the following. 
The float WL ambiguity can be derived from 

fiba — Pr; _ AP ra + foPro AT 


alias os O hèh ntn a 


where A; wL = c is the wavelength of WL ambiguity. The float WL ambiguity can be 
1 72 


further decomposed as 


—S 


Newt = New. + drwi — divı (10) 


fiDra + f2Dr,2 
Ar wi + fa) 


fi; Dj + fD? 
Ar. wis + f2) 


(11) 
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HMW combinations eliminate the geometric distance, satellite and receiver clock errors, 
as well as atmospheric delays. Averaging the HMW for a continuous arc results in accu- 
rate estimations of WL ambiguities, which can be easily fixed to integers by rounding. 


When the WL ambiguity Nw; is correctly resolved to an integer value and the float IF 


ambiguity N;. ir is obtained from PPP solution, the float NL ambiguity observable N,N L 
can be derived based on (8) 


—s fith, =s hr = 


r,NL — — Ae Ny ir = h-E rWL 7 Nr + dy NL vx dy; (13) 
where 
h tf 
d, NL = Tu (Ag By, ip = Dy tr) (14) 
ht fe 
dn, = T (Arr Bfr mi Dir) (15) 


Equations (10) and (13) serve as the basic model for estimating FCBs. Since they have the 
same structure, a general expression can be formulated as 


—S 


RÈ = N, — Ñ; = d, — d° (16) 
for WL and NL linear combinations, respectively. R7 represents the FCB measurements; 
N, denotes the float undifferenced ambiguities; ÑS denotes the integer part of N,, which 
is the sum of the original integer ambiguity and the integer part of the combined code 
hardware delays and uncalibrated phase delays from both receiver r and satellite s; d, 
and d denote the receiver and satellite FCBs. Note that N} can be determined by round- 
ing N? assuming the float ambiguities N? are precisely estimated. 


A set of equations in the form of (16) can be generated based on a network of reference 
stations. Ge et al. (2008) proposed to resolve the system by averaging based on all in- 
volved float single-differenced WL and NL ambiguity estimates. The single-differencing 
across satellites within single station eliminates the receiver FCB and single-differenced 
satellite FCB can be obtained. Since the same single-differencing approach will be imple- 
mented at the user end when using the FCB product, it is not necessary to obtain the 
absolute FCB and the single-differenced FCB would be sufficient. Based on the same 
principle but instead of the averaging process, a least squares method in an integrated 
adjustment was adopted to enhance the estimates (Li and Zhang 2012; Zhang and Li 
2013). Suppose that there are n satellites tracked at m reference stations, the system of 
equations can be expressed as 
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The obtained system of equations is singular. One arbitrarily selected FCB should be set 
to zero. In this way, the system of equations can be resolved with least squares algorithm 
(Li et al. 2015a). Note that the above discussions are based on single epoch data. Consid- 
ering the temporal stabilities of code hardware delays and uncalibrated carrier phase 
delays, the satellite FCBs can be estimated as constant values for specific time intervals, 
e.g., one day for WL FCB and 15 minutes for NL FCB (Ge et al. 2008). This can be accom- 
plished by accumulating all the equations within the specific time intervals. In this way, 
the final FCB products are produced and broadcasted to users. Since January 1, 2015, the 
School of Geodesy and Geomatics at Wuhan University (WHU-SGG) routinely releases 
GPS WL and NL FCB products with open access (Li et al. 2015a). 


With the above-estimated FCBs, users are able to correct the satellite FCBs, and the sin- 
gle-differenced PPP ambiguities are aimed to be fixed 


NP" = N,—N, — sm (18) 


For WL float ambiguities, they can be directly fixed by the rounding approach (Ge et al. 
2008). For NL float ambiguities, they are fed into the Least-squares AMBiguity 
Decorrelation Adjustment (LAMBDA) algorithm (Teunissen et al. 1997). If not all the 
float NL ambiguities can be fixed by the LAMBDA method, partial ambiguity resolution 
can be employed (Li and Zhang 2015; Teunissen et al. 1999). Once the WL and NL integer 
ambiguities are resolved and validated, a tight constraint can be reconstructed and im- 
posed on the ionospheric-free float ambiguities 


0=N,-N, -NS™—d5™, g, (19) 
In this way, a fixed solution can be obtained at the user end. 


The above discussed model is referred as FCB based PPP AR model. Collins et al. (2008) 
developed a decoupled clock model by separating satellite clocks for code and phase 
observations. Similarly, Laurichesse et al. (2009) developed an integer phase clock model 
in which the NL FCBs were assimilated into receiver and satellite clock estimates of a 
network solution. This model has been employed to generate the precise satellite clock 
products by Groupe de Recherche de Géodésie Spatiale of the Centre National d'Etudes 
Spatiales (CNES-GRG) (Loyer et al. 2012). The only difference between the decoupled 
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clock model and the integer phase clock model is the approach for determining the WL 
ambiguity. The integer phase clock model utilizes WL FCB corrections and a satellite- 
averaging process to fix the integer WL ambiguity, whereas the decoupled clock model 
directly estimates the integer WL ambiguity along with other unknowns through the 
functional model. The difference between the WL/NL FCB model and the integer phase 
clock model is the strategy of separating NL FCBs from integer ambiguities. NL FCBs 
are directly estimated in the WL/NL FCB model, whereas they are assimilated into the 
clock estimates in the integer phase clock model. These PPP AR techniques are proven 
to be equivalent in theory (Shi and Gao 2014). The systematic biases between position 
estimates have been demonstrated to be minimal and negligible with data from a global 
network of reference stations (Geng et al. 2010). The FCB based model can conveniently 
supplement current network solutions as an additional software module, as the FCB de- 
termination is compatible with current official clock-generation methods. For this rea- 
son, the FCB based model is adopted throughout the thesis. 


3 Problem statements and research objectives 


The PPP AR model and method described above have been proposed for decades. Due 
to the rapid development, and huge benefit, of multi-frequency and multi-GNSS meas- 
urements, it is required to extend the current model and method to integrate these meas- 
urements for fast and high precision positioning. 


3.1 Problem statements 


(1) Fast FCB estimation 


Based on the discussion in Section 2, we find that the least squares method is routinely 
utilized for FCB estimation (Li et al. 2015a; Li et al. 2017b; Liu et al. 2017; Wang et al. 
2017). However, the least squares method can be extremely time-consuming concerning 
the large number of observations from hundreds of reference stations and thousands of 
epochs during the FCB estimation. The dimension of the involved matrix could be fur- 
ther enlarged by multi-frequency and multi-GNSS observations, which also requires a 
lot of computation resource. This is unfavorable for the more and more popular real- 
time applications (Geng et al. 2018; Li et al. 2015b; Li et al. 2013c). In addition, iterations 
are required to deal with the one cycle inconsistency among FCB measurements. Since 
the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for WHU- 
SGG and CNES-GRG products respectively, the one cycle inconsistency arises whenever 
the superposition of receiver FCB and satellite FCB exceeds the boundary (Xiao et al. 
2018b). The additional iterations of the least squares method and computation of a large 
matrix demand a long time to finish. Therefore, a fast and efficient estimation method is 
desirable. 


Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution 


(2) Extend PPP AR to multi-GNSS 


Compared with GPS, the study of other GNSS PPP AR is limited and requires more at- 
tentions considering the rapid developments in recent years. For GLONASS, special care 
has to be taken to deal with the receiver inter-frequency bias between satellites originat- 
ing from the frequency division multiple access strategies. Usually, it is achieved by us- 
ing a network of homogeneous receivers (Geng and Bock 2016; Geng and Shi 2016; Liu 
et al. 2017). Also, BDS-2 code measurements suffer from the satellite induced pseudor- 
ange variations which are elevation- and frequency-dependent. Correction models based 
on multipath combination have been proven effective (Wanninger and Beer 2015). Un- 
fortunately, this approach does not apply to geostationary orbit satellites due to the al- 
most constant satellite elevation angle. Currently, only regional PPP AR or PPP AR with 
inclined geosynchronous orbit and medium earth orbit satellites is achievable for BDS 
(Gu et al. 2015a; Li et al. 2017a; Li et al. 2017b). Tegedor et al. (2016) initially assessed the 
Galileo PPP AR with four In-Orbit Validation satellites using regional stations around 
Europe. Li et al. (2017b) estimate the FCBs for the four global systems and assess the 
performance of combined PPP AR. The studies of Galileo and BDS PPP AR have been 
limited by the number of available satellites. The situation has been changed, with a lot 
of new Galileo and BDS satellites in the last two years (Xiao et al. 2019b). It is of interest 
to extend PPP AR to multi-GNSS, especially with the relatively new Galileo and BDS. 


(3) Extend PPP AR to multi-frequency 


The traditional PPP based on dual-frequency ionospheric-free combination cannot 
properly handle the additional signals, especially the functional and integer ambiguity 
resolution model. As a result, the PPP model based on uncombined measurements 
(Schónemann et al. 2011; Zhang et al. 2012), in which the individual signal of each carrier 
frequency is treated as an independent observable, has drawn increasing interest in the 
GNSS community (Gu et al. 2015b). Its efficiency has already been confirmed in terms of 
convergence time and precision for single-frequency PPP (Lou et al. 2016), multi-GNSS 
PPP (Chen et al. 2015), as well as PPP-RTK (Feng et al. 2013; Odijk et al. 2016). Moreover, 
this approach has been tested effectively for ionospheric modelling (Tu et al. 2013; Zhang 
et al. 2012), differential code bias (DCB) estimation (Liu et al. 2018; Shi et al. 2015), and 
satellite orbit determination (Strasser et al. 2018; Zehentner and Mayer-Gürr 2016). Com- 
pared to the well-developed IF PPP model, the uncombined PPP model, called uncom- 
bined PPP in the sequel, requires more investigations, especially in the case of multi- 
frequency processing and ambiguity resolution. First, how to model and constrain the 
ionospheric delay has a crucial impact on performance using the uncombined PPP. For 
example, it has been demonstrated that a white noise model is not adequate to capture 
the characteristics of the ionospheric delay. The external constraints developed from the 
ionospheric products, such as the IGS global ionosphere maps, are also not accurate 
enough to completely separate the ionospheric effects from the ambiguity parameters 
(Gu et al. 2015b). The influence of the ionospheric effects on the ambiguity fixing there- 
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fore must be reduced. Second, the method to deal with the DCB errors is more problem- 
atic with the uncombined PPP (Guo et al. 2015) than the dual-frequency IF PPP, since 
the latter can cancel out the DCB biases (Kouba and Héroux 2001). The problem of partial 
assimilation of the code bias (DCB) into phase bias (FCB) should also be carefully consid- 
ered. Third, the uncombined PPP approach was proposed to deal with multi-GNSS and 
multi-frequency signals, so a generalized FCB estimation and AR method (Li et al. 2018), 
which is extendable to dual-, triple-, and multi- frequency, should be proposed. 


Li et al. (2013b) verified the feasibility of the uncombined PPP AR with refined iono- 
spheric models. The ionospheric delay was constrained from a priori spatial-temporal 
information and ionospheric products. The GPS dual-frequency ambiguities were fixed 
sequentially in the forms of WL/NL, which followed the convention of IF PPP AR. Gu et 
al. (2015a) further testified the uncombined PPP AR with BDS triple-frequency observa- 
tions. The extra-wide-lane (EWL) and WL ambiguities were successfully fixed, whereas 
the B1 ambiguities were kept as float values. In addition, the performance was further 
limited by the satellite-induced multipath effects (Wanninger and Beer 2015). Li et al. 
(2018) proposed a unified FCB estimation and PPP AR method, which is extendable to 
multi-frequency uncombined PPP. The FCBs on each frequency were directly estimated 
from the raw float ambiguities derived from triple frequency observables. The model 
showed a great potential for multi-frequency uncombined PPP AR, although its DCB 
strategy may not be optimal (Xiao et al. 2019a). The satellite DCBs, together with the 
receiver DCB, were estimated as unknowns, and as a result the number of unknown 
parameters was increased. Given that the satellite GNSS DCB product is currently avail- 
able on a routine basis (Wang et al. 2016), it would be beneficial to make use of these 
products. 


In conclusion, it is required to develop a unified modeling strategy for multi-frequency 
uncombined PPP ambiguity resolution. 


(4) Data preprocessing in uncombined PPP 


The uncombined signal processing brings new challenges for cycle slip detection. The 
most significant feature is carrier frequency identification of cycle slips. Since all carrier 
frequency observations are processed separately, it is essential to identify the carrier fre- 
quency of the cycle slips to avoid contaminating other observations (Xiao et al. 2018a). 
In conventional methods, such as TurboEdit algorithm (Cai et al. 2013; Hatch 1982; Liu 
2011; Melbourne 1985; Wübbena 1985) and combinations of multi-frequency observa- 
tions (Dai et al. 2009; de Lacy et al. 2012; Wu et al. 2010a; Zhao et al. 2015), a common 
characteristic is the forming of optimal combinations to mitigate the presence of geomet- 
ric and ionospheric errors, which makes it impossible to identify the carrier frequency 
affected by a cycle slip. All carrier frequencies are flagged as containing cycle slips even 
if there is only one carrier frequency suffering from a cycle slip. Therefore, a cycle slip 
detection method, which is capable of carrier frequency identification of cycle slip and 
applicable to all GNSS should be proposed. 
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3.2 Research objectives 


The general research objective of this thesis is to improve the performance of GPS-only 
PPP AR by multi-GNSS integration and multi-frequency observations. In detail, the fol- 
lowing questions and aspects have been addressed within the thesis: 


e ls it possible to achieve PPP AR with the relatively new GNSS systems, e.g., Galileo 
and BDS? 


e How is the performance of PPP AR with the current constellation of Galileo and BDS? 
e How does PPP AR benefit from multi-GNSS integration? 


e Isit possible, and how, to achieve PPP AR with the multi-frequency observations, 
e.g., from Galileo and BDS? 


e How does PPP AR benefit from multi-frequency observations? 


e What is characteristic of FCB and how to efficiently estimate the FCB for real time 
PPP AR? 


e What is the uncombined PPP and what is the difference with regard to IF PPP? 
What influence can be expected, e.g., on data preprocessing, bias correction and 
parameter setting? 


4 Contributions 


4.1 Fast FCB estimation for PPP AR 


Estimating satellite FCBs based on the Kalman filter is proposed and demonstrated. 
Since the Kalman filter is based upon least squares methods (LSM), it is theoretically 
possible to calculate the same solution as for the commonly used LSM. In the proposed 
Kalman filter, the large number of observations is handled epoch-by-epoch, which sig- 
nificantly reduces the dimension of the involved matrix and accelerates the computation. 
Hence it outperforms the commonly used LSM in terms of efficiency. In order to avoid 
iterations caused by the one cycle inconsistency among FCB measurements, a pre-elim- 
ination method is developed based on the temporal stability of satellite FCBs. The pre- 
elimination method shows a clear advantage over post-residual adjustment, which fur- 
ther improves the efficiency (Xiao et al. 2018b). 


(1) Pre-elimination of one cycle inconsistency 


Since the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for 
WHU-SGG and CNES-GRGS products respectively, the one cycle inconsistency arises 
whenever the superposition of receiver FCB and satellite FCB exceeds the boundary. A 
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simple example is presented to depict the situation. Assuming there is one satellite (d? = 
0.2) tracked at two stations (d, = 0.6 and d, 2 = 0.8), the superposition of FCBs should 
be (Rz, = —0.4 and R}, = —0.6). However, the R7; would be 0.4 due to the rounding 
approach. As a consequence, the corresponding satellite FCB derived from the two sta- 
tions differs by one cycle. The one cycle inconsistency can be detected by examining the 
posterior residuals, as employed in former works. The posterior adjustments are ineffi- 
cient as iterations of the whole process are required. In contrast, we propose to eliminate 
the inconsistency in advance (shown in Figure 1), which shows a clear advantage over 
posterior adjustment. The proposed method consists of the following steps: 


1. For all the FCB measurements at a single station, the satellite FCBs are eliminated 
using previous estimates. The underlying assumption is that satellite FCBs are sta- 
ble over successive epochs and can be eliminated to a large extent by previous esti- 
mates. 

2. The residual of each FCB measurement after the subtraction of the satellite FCB 
yields an initial estimate of the receiver FCB. A set of initial receiver FCB estimates 
could be obtained from all simultaneously observed satellites. 

3. Intheory, the receiver FCBs obtained from all simultaneously observed satellites are 
expected to be consistent. Therefore, a one cycle inconsistency can be detected by 
examining the group of initial receiver FCBs. If an inconsistency exists, the corre- 
sponding receiver FCB will differ by +1 cycle. 

4. To compensate the one cycle inconsistency in receiver FCBs, the corresponding in- 
teger ambiguity N} is adjusted by +1 cycle. 

5. The procedures described above can be performed in an iterative way until a con- 
sistent receiver FCB is obtained from all satellite measurements. 


In addition, stations can also be handled individually to obtain clean FCB measurements. 
After the pre-elimination of the one cycle inconsistency, the clean FCB measurements 
are fed into the Kalman filter. 


(2) Design of Kalman filter for FCB estimation 


The Kalman filter addresses the general problem of state estimates of discrete time con- 
trolled processes that are governed by a linear stochastic difference equation (Kalman 
1960). The theory has been well studied and widely applied (Yang 2006; Yang et al. 2010). 
Since the Kalman filter is based upon the theory of least squares, it is theoretically pos- 
sible to calculate the same solution as LSM and versa vice. However, the design and 
normal equation matrix will be huge in LSM, considering the large number of observa- 
tions from hundreds of reference stations and thousands of epochs. Computation of the 
large matrix is time-consuming, which is inefficient and unfavorable for real-time appli- 
cations. In contrast, the large number of observations are handled epoch-by-epoch in a 
Kalman filter, which significantly reduces the dimension of the involved matrix and ac- 
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celerates the computation. One additional advantage of Kalman filter is its real-time ca- 
pability. This is of particular interest for real-time PPP AR and its applications. There- 
fore, a Kalman filter is employed in this work. The design and flowchart of a Kalman 
filter are depicted in Figure 2. In our case, R7 is the input FCB measurement for the Kal- 
man filter, while d, and d? are the output estimates. 


Receiver NL FCB [Cycle] Receiver WL FCB [Cycle] 


0 6 12 18 24 0 6 12 18 24 
UTC Time [Hour] UTC Time [Hour] 


Figure 1 Pre-elimination of the one cycle inconsistency for station IOOE on 2016/02/26. The top panels rep- 
resent receiver WL FCB adjustment, while the bottom ones represent that of receiver NL FCB. The num- 
bers denote GPS satellite PRNs. The left and right panels represent the raw and adjusted measurements, 
respectively. 
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Figure 2 Flow chart of the proposed Kalman filter. The dashed blocks represent additional steps adopted 
in this thesis. 


Two additional steps are added to the standard Kalman filter (Xiao et al. 2018b). The first 
step aims at establishing the dynamic model and determining the state transition matrix 
by analyzing the temporal stability of existing FCB products. Satellite WL FCBs are sta- 
ble over several days, which can be characterized as a constant parameter on a daily 
basis. Satellite NL FCBs are considered as constant within 15 minutes but exhibit small 
variations over 15-minute intervals. Therefore, NL FCBs are characterized as random 
walk processes. The second step is introduced to eliminate the one cycle inconsistency 
of FCB measurements. Only clean measurements are sent to the Kalman filter to avoid 
iterations. Note that a constraint is imposed on the Kalman filter to eliminate the rank 
deficiency. A satellite FCB is selected and set to zero, which is accomplished by adding 
a pseudo observational error equation with zero variance (Yang et al. 2010). 


4.2 Multi-GNSS PPP AR 


In the context of multi-GNSS, the PPP ambiguity resolution is extended to Galileo and 
BDS. In the view of the rapid developments of these systems in recent years, the feasi- 
bility of Galileo and BDS PPP ambiguity resolution with the current constellation is 
demonstrated. The satellite WL/NL FCBs are estimated from globally distributed Multi- 
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GNSS EXperiment (MGEX) stations and assessed by a comparison with those of GPS 
(Xiao et al. 2019b). 


(1) WL FCB comparison 


Results of 60 days indicate that the quality of Galileo WL FCB is better than for GPS and 
BDS in terms of data usage rate, residual distribution, as well as standard deviation 
(STD) of daily estimates (shown in Figure 3). We attribute the good quality of Galileo 
WL FCB to its advanced signal modulation, AltBOC, which significantly compresses the 
multipath effect for code measurements as shown Figure 4. 


0.04 - | HJ Foc | 
BEE lov BDS > 
0.03 | EEE Ecc L 
E! 
Ó GPS 
n 0.02 4 ^a [E 
> 
[65] 
NI | | | 
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PRN 


Figure 3 STDs of daily satellite WL FCBs with 60 days. Averaged STDs are calculated for GPS and BDS, 
while an individual STD is calculated for each Galileo satellite. E30 is selected as the reference satellite for 


Galileo, while G01 and C09 are selected for GPS and BDS, separately. 


1 T 
+ GAL 


MW [Cycle] 


20 40 60 80 
Elevation [Degree] 


Figure 4 HMW linear combinations of GPS, Galileo, and BDS dual frequency code measurements at sta- 


tion YEL2 on DoY 160, 2017. 
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Within the Galileo constellation, the quality of Full Operational Capability (FOC) satel- 
lite WL FCB is much better than for In-Orbit Validation (IOV) satellites. The root mean 
square (RMS) for each satellite in Galileo are presented in Figure 5. It is found that the 
RMS of IOV satellites are significantly larger than other FOC satellites, while there is no 
significant difference between FOC and two highly eccentric (ECC) satellites. These re- 
sults indicate that the quality of FOC WL float ambiguities is better than that of IOV 
satellites. The worse performance of IOV satellites probably stems from the reduction of 
signal transmission power for IOV satellites. ESA imposed a reduction of 1.5 dB in the 
signal power of all four IOV satellites following a payload power problem of the fourth 
IOV (E20) in 2014. The 1.5 dB power decrease results in an increase of approximately 15- 
20% of the thermal noise of the receiver, which roughly matches the observed increase 
in IOV WL FCB residuals. 
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Figure 5 RMS of WL residuals after WL FCB estimations. Averaged RMSs are calculated for GPS and BDS, 
while an individual RMS is calculated for each Galileo satellite. 


(2) NL FCB comparison 


Figure 6 depicts the usage rate of float ambiguities for each Galileo satellite (average 
94.0%), as well as the average usage rates for GPS (97.6%) and BDS (80.4%). The usage 
rate of GPS NL float ambiguities is highest while that of BDS is worst. Since the weights 
of carrier phase measurements are 10000 times larger than those of code measurements, 
the quality of the float ambiguities is dominated by the unmodeled errors in PPP. There- 
fore, the highest usage of GPS is reasonable since its precise product and models are 
currently the state of the art. All Galileo satellites except E19 show similar usage rates 
ranging from 91.8 to 96.9%, indicating a similar quality level of NL ambiguities. The us- 
age rate of E19 is only 83.876, which is much smaller than that for other Galileo satellites. 
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Inspired by Zaminpardaz and Teunissen (2017) who reported that for SEPTENTRIO re- 
ceivers the carrier-to-noise density ratio for E19 lies below the value of the other two IOV 
satellites for elevations higher than 60 degrees, we further calculated the carrier-to-noise 
density ratio, data usage rates for three receiver groups, namely LEICA (37 stations), 
SEPTENTRIO (41 stations) and TRIMBLE (83 stations). From the results in Figure 7, it 
can be seen that the carrier-to-noise density ratio of E19 lies below the other satellites for 
all the receiver groups. In addition, the data usage rates of E19 are lower than for the 
other satellites for all the three groups and no clear relationship with receiver type is 
found. Furthermore, we notice that the STD of satellite laser ranging (SLR) residuals for 
E19 (0.130 m) is larger than for E11 (0.092 m) and E12 (0.086 m) when performing SLR 
validation (Guo et al. 2017). This may indicate a poor quality of precise satellite products 
for E19. Based on the above analysis, we suspect that the lower data usage rate of E19 is 
caused by the joint effect of poorer satellite products and larger carrier phase noise. 
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Figure 6 Usage rates of float NL ambiguities for satellite NL FCB estimation. Averaged usage rates are cal- 
culated for GPS and BDS, while an individual usage rate is calculated for each Galileo satellite. 
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Figure 7 Carrier-to-noise density ratio C/NO of Galileo signals for three receiver groups, TRIMBLE (left), 
SEPTENTRIO (middle) and LEICA (right). 
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In conclusion, the quality of Galileo NL FCB is slightly worse than that of GPS but better 
than that of BDS. Within Galileo, the NL FCB quality of FOC and IOV (except E19), as 
well as the two eccentric satellites, shows no significant difference in terms of data usage 
rates and residuals. The reason for the worse performance of E19 is still not clear. On the 
one hand, it cannot, or at least not fully, be ascribed to the signal transmission power as 
the power difference between FOC and IOV E11/E12 is larger than the difference be- 
tween E11/E12 and E19. On the other hand, the worse quality of E19 satellite orbit, indi- 
cated by SLR residuals, could also degrade NL FCB estimations. This issue remains an 
open question and deserves further investigation. 


For static Galileo PPP AR with three-hour sessions, the positional biases can be reduced 
by 67, 45 and 22% for east, north and up components, respectively. In addition, PPP AR 
also improves the kinematic solutions by 59, 45 and 28% for east, north and vertical com- 
ponents, respectively. 


4.3 Multi-frequency PPP AR 


A unified model for multi-frequency PPP AR based on raw uncombined observations is 
proposed (Xiao et al. 2019a), which simplifies the concept of phase biases for AR. No 
assumption is made on the method used to determine FCB on the server end, which 
implies that the generated FCB from multi-frequency observations could be flexibly 
used, such as for dual- or triple-frequency ones. It is demonstrated that the model is 
extendable to dual- and triple-frequency observations. 


(1) Multi-frequency uncombined PPP AR model 


To deal with multi-frequency observations, the uncombined PPP model is adopted. In 
the uncombined PPP model, the ionospheric delay is directly estimated. Another im- 
portant difference between the IF PPP and the uncombined PPP models is the strategy 
to deal with the DCB. The DCB is not of concern in IF PPP as the IF combination is also 
used for precise clock generation, which implies that the DCB could be fully absorbed 
by other parameters (Zhang et al. 2012) or simply cancelled out in the IF PPP (Dach et al. 
2015). But this is not the case in the uncombined PPP, especially with multi-frequency 
observations. We propose estimation of the receiver DCB and correction of the satellite 
DCB with existing multi-GNSS DCB products (Wang et al. 2016). Taking triple-frequency 
observations as an example, the correction equation for the raw pseudorange observa- 
tions can be deduced as (Guo et al. 2015) 


= 1 
27 


- a 
Pro = pp + dt, — dtre + dT + az ` dlp, — DCBi + DCB,,5 + £p, (20) 
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_ 1 
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where DCB§, = Dj — Di, DCB,42 = Dy; — D,4, and dt, = dt, + D,,. The DCB$, and 
DCB;3 can be obtained from multi-GNSS DCB products, while DCB, 12 and DCB, 13 are 
estimated as daily constant parameters. Similarly, the phase equations can be rewritten 
as 


np = Ppt dt,— dtju, + dT — aj: dl2, + AN}, + Eo, (21) 


where the ambiguity can be re-parameterized as 
775 E S 
Ny f = Nèf + b, f = be 
byf = B, f z D,4A/Ay (22) 


and the estimable parameters are 
-- s —S —$ —s 
X= |x y z dt, dT Try Dra Dri3 Nra Ny2 LAM (23) 


The estimated ambiguity parameter is a combination of the integer ambiguity, the cor- 
responding code hardware delays, and the uncalibrated carrier phase delays at both re- 
ceiver and satellite ends. Similar to that of IF PPP, the integer property of the ambiguity 
parameter can be recovered provided FCB is corrected. In dual-frequency IF PPP, the 
float ambiguity is usually decomposed into WL/NL forms in order to recover the integer 
property (Ge et al. 2008). This is partly because the IF combination of L1/L2 ambiguities 
is, in essence, not an integer. Another reason is that the WL ambiguities possess a rela- 
tively longer wavelength and are less correlated, therefore can be easily fixed. For un- 
combined PPP AR, it is also important to form combinations of raw ambiguities. On the 
one hand, the estimated raw float ambiguities are strongly correlated (Li et al. 2018). On 
the other hand, the raw float ambiguities are quite sensitive to unmodeled ionospheric 
errors (Gu et al. 2015b). Therefore, the combinations with longer wavelengths and lower 
ionospheric delays are preferred. From the systematic investigation of triple-frequency 
combinations, it is found that the following combinations are a good compromise be- 
tween ionospheric reduction and noise amplification. These combinations are used for 
both BDS and Galileo triple-frequency observations for simplicity 


—S —S 


Nri 4 —3 0 Nra 

Nr ic2| = h =l ] N,2 (24) 
—S 1 0 -11|—s 

N. 1c3 Nx3 


Substituting (22) into the above system produces the basic model for estimating FCBs. 
Now the equations are quite similar to (16), which implies that the same method can be 
used to estimate FCB for these combinations, i.e., (17). 
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With the obtained combined FCB, we are able to calculate the FCB of the raw L1/L2/L3 
carrier frequency 


d| [4 -3 oOrp'![dic 
d$ -h =j ] di (25) 
dj l 0 -1l [des 


The transformation from combined FCBs to raw FCBs is important, as it provides more 
flexibility to the users. With the raw FCBs, users are able to choose their own linear com- 
binations of observations, formulate the corresponding combined FCB, and conduct PPP 
AR. This representation allows interoperability if the server and user sides implement 
different AR methods. In addition, the raw FCB is suitable for the State Space Represen- 
tation (SSR) of Radio Technical Commission for Maritime services (RTCM) (Weber et al. 
2005), where one phase bias per phase observable is broadcasted instead of making spe- 
cific combinations. 


Similar to dual-frequency IF PPP AR, single differencing across satellites must be firstly 
performed in order to remove receiver FCBs. Then the single-differenced ambiguities 
from different carrier frequencies are combined, as has been done during FCB estimation 


mn mn 
Nia u A dk] 

j à r apn 
N, ijk2 i; j2 k2||N,2 (26) 
arm,n i 1 k am,n 
N, ijk3 3 Js Rs] Ns 


where N"? = N. — N. is the single-differenced ambiguity between satellites m and n. 
Based on the coefficients (i, j,k), the FCB for the specific combined ambiguity can also 
be formed. Note that the linear combinations are not necessarily to be the same as those 
in FCB estimation, although the three combinations mentioned above are strongly rec- 
ommended. In our experiments, we have used the same combinations as in FCB gener- 
ation for uncombined PPP AR. Usually, the EWL/WL float ambiguities can be directly 
fixed by the rounding approach after the correction of FCB (Ge et al. 2008), and the NL 
float ambiguities are fed into the LAMBDA algorithm to search for correct integers 
(Teunissen et al. 1997). However, in our study, the LAMBDA is used for each combina- 
tion, regardless of its property, which simplifies the design of the algorithm. In addition, 
if not all the float ambiguities can be fixed by the LAMBDA method, partial ambiguity 
resolution can be employed (Li and Zhang 2015; Teunissen et al. 1999). It is found that 
the searching and fixing of ambiguities for the combination with longer wavelengths 
(e.g., EWL/WL) is quite fast. When the integer ambiguities for one combination are re- 
solved and validated, a tight constraint can be reconstructed. The number of constraints 
accumulate as the process repeats for all linear combinations. Afterwards, the constraints 
are imposed on the raw ambiguities, and yields the AR solution. Note that the ambigui- 
ties in IF PPP AR must be sequentially fixed in the order of WL/NL. An IF ambiguity is 
constrained only when both its WL and NL ambiguities are fixed, while the linear com- 
bined ambiguities in our approach can be fixed and constrained independently. 
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(2) Test results with Galileo and BDS triple-frequency observations 


To verify the efficiency of the proposed approach, we processed 51 days of Galileo and 
BDS triple-frequency observations collected from globally distributed MGEX stations 
(Xiao et al. 2019a). The estimated FCB shows a good consistency with the input float 
ambiguities. The RMS of Galileo FCB residuals is 0.05 cycles, while that of BDS is 0.08 
cycles. It is also observed that the residuals are smaller for the combinations with larger 
wavelengths. The results indicate that there may exist ionospheric errors, and combina- 
tions are required to reduce their influence. The average STD of combined FCB is around 
0.03 cycles, while that for raw FCB is around 0.10 cycles, as shown in Figure 8 and Figure 
9. To reduce the communication with servers for real time applications, it would be more 
efficient to broadcast linear combined FCB. When comparing the results from multi- 
GNSS, it can be seen that the STDs of Galileo EWL/WL FCBs are smaller than those of 
GPS and BDS, while that of the Galileo NL FCB is worse than GPS and BDS. The better 
quality of Galileo EWL/WL FCBs is likely attributed to the multipath suppression of Gal- 
ileo signals, while the worse quality of Galileo NL FCB is due to the poor precision of 
satellite orbit and clock product. It is found that the STDs of Galileo raw FCBs are smaller 
than that of GPS and BDS, regardless of combinations for most of the days. 


Mean STD [cycle] 


Figure 8 Mean STD of the combined FCB series for all 51 days. Daily STD is calculated for each satellite 
FCB series. For each day, the mean STD of all satellite daily STDs is presented. 
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Mean STD [cycle] 


Figure 9 Mean STD of the raw FCB series for all 51 days. Daily STD is calculated for each satellite FCB se- 
ries. For each day, the mean STD of all satellite daily STDs is presented. 


The performance of triple-frequency PPP AR is assessed with 11 days of data in three- 
hour sessions. Compared with the float solutions, the positional biases of AR solutions 
are significantly improved, as shown in Figure 10 and Figure 11. The improvements of 
east, north, up components for positioning with BDS are 49.2%, 38.3%, and 29.6%, while 
those for positioning with Galileo are 60.0%, 29.0%, and 21.1%. These results demon- 
strate the efficiency of the proposed FCB estimation approach, and that the triple-fre- 
quency PPP AR can bring an obvious benefit to the float solution. 
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Figure 10 Convergence performance of BDS triple-frequency PPP float and AR solutions based on 804 3-h 
sessions under 68% confidence level. 
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Figure 11 Convergence performance of Galileo triple-frequency PPP float and AR solutions based on 5805 
3-h sessions under 68% confidence level. 


When comparing triple-frequency PPP AR with that of dual-frequency, it is found that 
the contribution of the third frequency observations is not remarkable. The insignificant 
improvement of the third frequency observable may be due to its limited contribution to 
satellite geometry and the narrow deployment with respect to the second carrier fre- 
quency. Nevertheless, adding the third frequency increases the reliability since it is ob- 
served that the number of successful sessions is increased. It is noted that all the experi- 
ments were conducted with excellent observational condition, e.g., the state-of-the-art 
receiver and antenna, open environment and static mode. Under this circumstance, the 
dual-frequency observations are of high quality and sufficient for high-precision posi- 
tioning, which also limits the contribution of the additional signals. More obvious im- 
provements could be expected for low cost receiver and antenna, restricted environment 
or kinematic mode (Li 2018). 


4.4 Uncombined cycle slip detection and repair 


As has been proved in the above research, the uncombined PPP is a good choice for the 
integration of multi-frequency and multi-GNSS observations. However, the uncombined 
signal processing brings new challenges for cycle slip detection, which is the carrier 
frequency identification of cycle slip (Xiao et al. 2018a). We presented an improved ap- 
proach based on a time-differenced model for cycle slip detection and repair in uncom- 
bined PPP. 


For cycle slip detection, the proposed approach significantly reduces false alarms as 
shown in Figure 12 and Figure 13. Having access to a reliable cycle slip detection method 
greatly reduces the number of ambiguity parameters to be estimated for processing un- 
combined GNSS measurements. 
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Figure 12 Results of cycle slip detection for rover (left) and reference stations (right) using geometry-free 
combinations (1, 1, -2) and (1, 2, 1) (Wu et al. 2010b; Zhen et al. 2008). Red triangle, green square, and blue 
inverted triangle represent B1, B2 and B3 BDS carrier frequency observations, respectively. 
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Figure 13 Results of cycle slip detection for rover (left) and reference (right) stations using the improved 
algorithm. 


Compared with the routine time-differenced methods (Banville and Langley 2013; 
Zhang and Li 2016), the proposed approach can not only identify the carrier frequency 
in which the cycle slip occurs, but also makes it possible to separate the observation at 
other frequency of the same satellite without cycle slip, which greatly contributes to cycle 
slip estimations. Simulation results show that the theoretical success rates increase to 
99.99% for both dual and triple frequency observations, as shown in Figure 14. Results 
from a real dataset also indicate that much stronger model strength of cycle slip estima- 
tion can be achieved. 
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Figure 14 The ambiguity dilution of precision (ADOP) based success rates of dual-frequency (left) and tri- 


ple-frequency (right) observations. The yellow surface represents the routine time-differenced model, 
while the green one represents the improved algorithm. 


Although the results shown above are obtained with BDS triple-frequency observations, 
the method has been implemented in our software and used for all the researches con- 


ducted in this thesis. It is proved that the method is easily applicable to dual-frequency 
and multi-frequency observations for all GNSS. 


4.5 Program development 


Program development is also an important contribution of this thesis. Two software 
suites, i.e. PPP ambiguity resolution engine and FCB estimation server, have been cre- 


ated. The PPP ambiguity resolution engine is developed based on RTKLIB 
(http://www.rtklib.com/), which is an open source program package for standard and 
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precise positioning with GNSS written in C programming language. The RTKLIB pack- 
age includes a portable program library and several application programs. The original 
RTKLIB provides excellent support for GNSS data processing, e.g., standard data format 
and basic functions for PPP. During my research, I have significantly extended the orig- 
inal RTKLIB software. These improvements mainly include: 


(1) Fixing bugs in the original software, for example, the stochastic model, antenna 
models and DCB corrections. 


(2) Enhancing the data pre-processing, e.g., cycle slip detection and repair, BDS-2 
satellite induced multipath correction. 


(3) Implementing FCB- and integer phase clock-based ambiguity resolution for 
dual-frequency ionospheric-free PPP, which includes reading FCB files, correct- 
ing satellite FCBs and fixing the ambiguity with the LAMBDA method. Partial 
ambiguity resolution may be employed when not all float ambiguities can be 
fixed. 


(4) Implementing the multi-frequency uncombined PPP model and ambiguity res- 
olution. For multi-frequency observations, the uncombined PPP model is imple- 
mented, in which the slant ionospheric delay and receiver DCB are directly esti- 
mated. Furthermore, combinations of multi-frequency raw ambiguities are 
formed and fixed. 


In addition, an FCB estimation software is developed with Python programming lan- 
guage, accompanied by a FCB analysis tool based on MATLAB programming language. 
The software supports FCB generation for both dual-frequency ionospheric-free PPP and 
multi-frequency uncombined PPP models. In order to investigate the characteristics of 
the FCB product, a detail report, which includes quality control information and residu- 
als statistics, is also generated. Besides the two software suites, various scripts have been 
developed for the efficient processing of huge GNSS data sets, such as automatic data 
preparing, batch processing and analysis of results. 


5 Publication overview 


This thesis includes the following four peer-reviewed publications (Xiao et al. 2019a; 
Xiao et al. 2019b; Xiao et al. 2018a; Xiao et al. 2018b): 


5) Xiao G, Li P, Gao Y, Heck B (2019a) A Unified Model for Multi-Frequency PPP 
Ambiguity Resolution and Test Results with Galileo and BeiDou Triple-Fre- 
quency Observations. Remote Sensing 11(2):116 doi:10.3390/rs11020116 


Author's contributions: The first author derived the unified model for multi-frequency 
PPP AR. The first and second author conceived and designed the experiments. The first 
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author collected the data, performed the experiments, analyzed the results, and wrote 
the paper. All authors reviewed the paper. 


6) Xiao G, Li P, Sui L, Heck B, Schuh H (2019b) Estimating and assessing Galileo 
satellite fractional cycle bias for PPP ambiguity resolution. GPS Solutions 23:3 
doi:10.1007/s10291-018-0793-z 


Author's contributions: The first author conceived the presented idea. The first and sec- 
ond author designed the experiments, worked out the technical details. The first author 
performed the numerical calculations for the experiments, analyzed the data, and wrote 
the paper. All authors reviewed the paper. 


7) Xiao G, Sui L, Heck B, Zeng T, Tian Y (2018b) Estimating satellite phase fractional 
cycle biases based on Kalman filter. GPS Solutions 22:82 doi:10.1007/s10291-018- 
0749-3 


Author's contributions: The first author devised the main conceptual ideas and de- 
signed the proposed Kalman filter for FCB estimation. The first author performed the 
experiments, analyzed the data, and wrote the paper. All authors reviewed the paper 
and helped shape the research. 


8) Xiao G, Mayer M, Heck B, Sui L, Zeng T, Zhao D (2018a) Improved time-differ- 
enced cycle slip detect and repair for GNSS undifferenced observations. GPS So- 
lutions 22:6 doi:10.1007/s10291-017-0677-7 


Author's contributions: The first author designed the model and the experiments. The 
first author performed the experiments, analyzed the data, and wrote the paper. All au- 
thors discussed the results, reviewed and contributed to the final manuscript. 


6 Conclusions and outlook 


6.1 Conclusions 


For decades, GPS PPP has found many scientific and industrial applications due to its 
cost-effectiveness, global coverage, and high accuracy. However, it suffers from a few 
drawbacks which prevents more applications, e.g., slow convergence, availability and 
reliability. The rapid development of multi-GNSS and multi-frequency signals provides 
an excellent opportunity for improvements. In addition, the integer ambiguity resolution 
technique is expected to further enhance the accuracy and shorten the convergence time. 
In the framework of this thesis, multi-frequency and multi-GNSS measurements are 
properly integrated for PPP with ambiguity resolution to achieve state-of-the-art fast and 
accurate positioning, which provides an important contribution to GNSS precise posi- 
tioning and applications. 
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Integer ambiguity resolution is the key issue to improve the PPP positioning quality. The 
core of PPP ambiguity resolution is carrier phase biases estimation. In this thesis, the 
characteristics of phase fractional cycle biases are analyzed with the updated and long 
time series. Results indicate that the stability of phase fractional cycle biases has been 
improved in recent years. Taking advantage of this property, an improved FCB estima- 
tion method based on Kalman filter is proposed. In the proposed Kalman filter, the large 
number of observations is handled epoch-by-epoch, which significantly reduces the di- 
mension of the involved matrix and accelerates the computation. Hence it outperforms 
the commonly used least squares method in terms of efficiency. In order to avoid itera- 
tions caused by the one cycle inconsistency among FCB measurements, a pre-elimination 
method is developed based on the temporal stabilities of satellite FCBs. The pre-elimi- 
nation method shows a clear advantage over post-residual adjustment, which further 
improves the efficiency. The fast estimation of FCB is especially useful for real-time ap- 
plication. Results from about 200 IGS stations indicate that the determined GPS satellite 
FCBs show a good consistency with existing FCB products, e.g., from IGS analysis center 
CNES-GRG and Wuhan University. 


With the estimated FCB products, the PPP ambiguity resolution technique with the cur- 
rent Galileo and BDS constellation is verified. The satellite fractional cycle biases are 
thoroughly assessed by a comparison with those of GPS in terms of data usage rate, re- 
sidual distribution, as well as standard deviation of daily estimates. Results indicate that 
the quality of Galileo WL FCB is better than for GPS and BDS. We attribute the good 
quality of Galileo WL FCB to its advanced signal modulation, AltBOC, which signifi- 
cantly compresses the multipath effect for code measurements. Within the Galileo con- 
stellation, the quality of FOC WL FCB is much better than for IOV satellites. The poorer 
performance of IOV satellites WL FCB is a result of a reduction in the satellite transmis- 
sion signal power. The performance of the two satellites with highly eccentric orbits is 
comparable to other FOC satellites but having a smaller number of observations. As for 
NL FCB, the quality of Galileo NL FCB is slightly worse than that of GPS but better than 
that of BDS. Since the accuracy of NL FCB estimation is dominated by unmodeled errors 
in float PPP, the worse quality of Galileo NL FCB is likely caused by the imperfect an- 
tenna models. Within Galileo, the NL FCB quality of FOC and IOV (except E19), as well 
as the two eccentric satellites, shows no significant difference in terms of data usage rates 
and residuals. The reason for the worse performance of E19 is still not clear. On the one 
hand, it cannot, or at least not fully, be ascribed to the signal transmission power as the 
power difference between FOC and IOV E11/E12 is larger than the difference between 
E11/E12 and E19. On the other hand, the worse quality of the E19 satellite orbit, indicated 
by SLR residuals, could also degrade NL FCB estimations. This issue remains an open 
question and deserves further investigation. Galileo PPP AR solutions are conducted at 
20 MGEX stations with three-hour sessions for ten days. The positional biases of AR so- 
lutions are 0.7, 0.6 and 2.1 cm for east, north and up components respectively, while 
those for float solutions are 2.1, 1.1 and 2.7 cm, corresponding to improvements of 67, 45 
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and 22% respectively. These results demonstrate that the Galileo observations can bring 
an obvious benefit to ambiguity-fixed PPP. 


PPP ambiguity resolution is also extended to multi-frequency observations in this thesis. 
In order to properly integrate the multi-frequency observations, a unified uncombined 
PPP ambiguity resolution model based on raw observations is proposed. To verify its 
efficiency, we processed 51 days of Galileo and BDS triple-frequency observations col- 
lected from globally distributed MGEX stations. The estimated FCB show a good con- 
sistency with the input float ambiguities. The RMS of Galileo FCB residuals is 0.05 cycles, 
while that of BDS is 0.08 cycles. It is also observed that the residuals are smaller for the 
combinations with larger wavelengths. The results indicate that there may exist iono- 
spheric errors, and linear combinations are required to reduce their influence. The aver- 
age STD of combined FCB is around 0.03 cycles, while that for raw FCB is around 0.10 
cycles. To reduce the communication with servers for real time applications, it would be 
more efficient to broadcast linear combined FCB. The performance of triple-frequency 
PPP AR is assessed with 11 days of data in three-hour sessions. Compared with the float 
solutions, the positional biases of AR solutions are significantly improved. The improve- 
ments of ENU components for positioning with BDS are 49.2%, 38.3%, and 29.6%, while 
those for positioning with Galileo are 60.0%, 29.0%, and 21.1%. These results demon- 
strate the efficiency of the proposed FCB estimation approach, and that the triple-fre- 
quency PPP AR can bring an obvious benefit to the float solution. When comparing tri- 
ple-frequency PPP AR with that of dual-frequency, it is found that the contribution of 
the third frequency observations is minimal. The insignificant improvement of the third 
frequency observable may be due to its limited contribution to satellite geometry and 
the narrow deployment with respect to the second carrier frequency. Nevertheless, add- 
ing the third frequency increases the reliability since it is observed that the number of 
successful sessions is increased. 


The uncombined signal processing brings new challenges for cycle slip detection. The 
most significant feature is carrier frequency identification of cycle slips. Since all carrier 
frequency observations are processed separately, it is essential to identify the carrier fre- 
quency of the cycle slips to avoid contaminating other observations. To provide contin- 
uous carrier phase measurements for the above research, an improved approach based 
on a time-differenced model for cycle slip detection and repair is proposed, which re- 
duces the false alarms and increases the success rate of cycle slip estimation. The simu- 
lation results show that the theoretical success rates increase to 99.99% for both dual and 
triple frequency observations. Results from a real dataset also indicate that much 
stronger model strength of cycle slip estimation can be achieved. 


In summary, the GPS only PPP AR has been extended to multi-frequency and multi- 
GNSS PPP ambiguity resolution with improved accuracy and fast convergence in this 
thesis. This is accomplished by a unified model based on the uncombined PPP. The pro- 
posed model has been carefully studied and enriched with improved cycle slip detection 
and repair, fast FCB estimation and ambiguity resolution methods. 
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6.2 


Chapter 1 Introductory chapter 


Outlook 


The outlook and perspectives of this thesis are presented as follows: 


1. 


Performance evaluation with full operational constellation. In the next years, 
Galileo and BDS will provide full operational service (Montenbruck et al. 2017). 
A significant increase of the number of usable satellites and signals can be ex- 
pected (Tian et al. 2019). It would be interesting to apply the model proposed in 
this thesis to the full operational GNSS and to evaluate the performance. 


Refinements of antenna corrections. Although the performance of Galileo PPP 
AR is still worse than that of GPS, a promising improvement can be expected in 
the near future as the Galileo FOC satellite metadata, including satellite mass, 
attitude law, PCO and PCV, are published. This information is expected to im- 
prove the precise orbit determination and PPP solution (Li et al. 2019). 


FCB broadcasting scheme. The traditional FCB scheme, which are broadcasted 
in WL and NL forms at 24-hour and 15-minute intervals (Ge et al. 2008) , respec- 
tively, needs to be updated. First, the stability has been improved with the ad- 
vancements of error models (Xiao et al. 2018b). Second, it is not suitable for 
multi-frequency uncombined PPP AR (Xiao et al. 2019a). Therefore, designing a 
novel scheme based on the characteristics of the uncombined model and the cur- 
rent FCB is demanding. 


Optimal combinations for multi-frequency ambiguity resolution. Tradition- 
ally, the IF ambiguities are decomposed into WL and NL ambiguities, which are 
sequentially fixed in IF PPP AR (Ge et al. 2008; Li et al. 2013b). However, in 
multi-frequency uncombined PPP AR, there are maybe more suitable combina- 
tions. The LAMBDA method can be used to automatically search for the optimal 
linear combinations of ambiguities (Li et al. 2018). With huge data sets, the opti- 
mal combinations can be found and recommended to users. 
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Abstract Phase fractional cycle biases (FCBs) originating from satellites and receivers destroy the 
integer nature of PPP carrier phase ambiguities. In order to achieve integer ambiguity resolution 
of PPP, FCBs of satellites are required. In former work, least squares methods (LSM) are com- 
monly adopted to isolate FCBs from a network of reference stations. However, it can be extremely 
time-consuming concerning the large number of observations from hundreds of stations and 
thousands of epochs. In addition, iterations are required to deal with the one cycle inconsistency 
among FCB measurements. We propose to estimate the FCB based on a Kalman filter. The large 
number of observations are handled epoch by epoch, which significantly reduces the dimension 
of the involved matrix and accelerates the computation. In addition, it is also suitable for real- 
time applications. As for the one cycle inconsistency, a pre-elimination method is developed to 
avoid iterations and posterior adjustments. A globally distributed network consisting of about 
200 IGS stations is selected to determine the GPS satellite FCBs. Observations recorded from DoY 
52 to 61 in 2016 are processed to verify the proposed approach. The RMS of wide lane (WL) pos- 
terior residuals is 0.09 cycles while that of the narrow lane (NL) is about 0.05 cycles, which indi- 
cates a good internal accuracy. The estimated WL FCBs also have a good consistency with existing 
WL FCB products (e.g., CNES-GRG, WHU-SGG). The RMS of differences with respect to GRG 
and SGG products are 0.03 and 0.05 cycles. For satellite NL FCB estimates, 97.9% of the differences 
with respect to SGG products are within + 0.1 cycles. The RMS of the difference is 0.05 cycles. 
These results prove the efficiency of the proposed approach. 


Keywords Precise point positioning : Integer ambiguity resolution - Fractional cycle bias - Kalman 
filter 


1 Introduction 


Integer ambiguity resolution (AR) can significantly shorten the convergence time and 
improve the accuracy of Precise Point Positioning (PPP). Phase fractional cycle biases 
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(FCBs) originating from satellites and receivers destroy the integer nature of PPP carrier 
phase ambiguities. The receiver FCB can be eliminated by single differencing across sat- 
ellites or assimilated into the receiver clock parameter by forcing one zero difference 
ambiguity to its nearest integer, while the satellite FCBs must be estimated from a net- 
work of reference stations (Gabor and Nerem 1999). With the additional satellite FCB 
estimates, PPP users are able to remove satellite FCBs and recover the integer nature of 
ambiguities. 


The fact that double-differenced ambiguities in global or regional networks can be re- 
solved to integer values lays the foundation of integer ambiguity resolution for PPP. The 
resolved double-differenced ambiguity implies that the fractional parts of two single- 
differenced ambiguities (across satellites) must agree well with each other. By estimating 
the fractional parts at the server end and applying them to single differencing PPP at the 
user end, PPP integer ambiguity resolution can be achieved. In this sense, PPP integer 
ambiguities are in fact equivalent to double-differenced ambiguities (Teunissen and 
Khodabandeh 2015). Ge et al. (2008) proposed to estimate the common fractional parts 
of the single-differenced ambiguities in wide lane (WL) and narrow lane (NL) form. The 
fractional parts, denoted as single-differenced FCB (SD-FCB), were estimated by averag- 
ing the fractional parts of all involved float single-differenced WL and NL ambiguity 
estimates. Based on the same principle but instead of the averaging process, a least 
squares method in an integrated adjustment was adopted to enhance the estimates (Li 
and Zhang 2012; Zhang and Li 2013). Since January 1, 2015, the School of Geodesy and 
Geomatics at Wuhan University (WHU-SGG) routinely releases GPS WL and NL FCB 
products with open access (Li et al. 2015a). Similar approaches have also been applied to 
BDS (Li et al. 2017a; Wang et al. 2017), Galileo (Tegedor et al. 2016) and GLONASS (Geng 
and Shi 2016). In order to exploit the ionosphere characteristics, the model has been ex- 
tended to deal with GPS L1 and L2 raw observables (Gu et al. 2015b; Li et al. 2013a) and 
BDS triple-frequency observables (Gu et al. 2015a). 


Different from the above approaches, Collins et al. (2008) developed a decoupled clock 
model by separating satellite clocks for code and phase observations. Similarly, Lau- 
richesse et al. (2009) developed an integer phase clock model in which the NL FCBs were 
assimilated into receiver and satellite clock estimates of a network solution. This model 
has been employed to generate the precise satellite clock products by Groupe de Recher- 
che de Géodésie Spatiale of the Centre National d'Etudes Spatiales (CNES-GRG) (Loyer 
et al. 2012). The only difference between the decoupled clock model and the integer 
phase clock model is the approach for determining the WL ambiguity. The integer phase 
clock model utilizes WL FCB corrections and a satellite-averaging process to fix the in- 
teger WL ambiguity, whereas the decoupled clock model directly estimates the integer 
WL ambiguity along with other unknowns through the functional model. The difference 
between the WL/NL FCB model and the integer phase clock model is the strategy of 
separating NL FCBs from integer ambiguities. NL FCBs are directly estimated in the 
WL/NL FCB model, whereas they are assimilated into the clock estimates in the integer 
phase clock model. 
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These PPP AR techniques are proven to be equivalent in theory (Shi and Gao 2014). The 
systematic biases between position estimates have been demonstrated to be minimal and 
negligible with data from a global network of reference stations (Geng et al. 2010). The 
WL/NL FCB model can conveniently supplement current network solutions as an addi- 
tional software module, as the FCB determination is compatible with current official 
clock-generation methods. In contrast, the integer phase clock products are incompatible 
with current clock products although they may perform slightly better in practice. 


Based on the review of existing work, we find that the least squares method (LSM) is 
routinely utilized for FCB estimation. However, the LSM can be extremely time-consum- 
ing concerning the large number of observations from hundreds of reference stations 
and thousands of epochs during the FCB estimation. This is unfavorable for the more 
and more popular real-time applications. In addition, iterations are required to deal with 
the one cycle inconsistency among FCB measurements. Since the FCB estimates are lim- 
ited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for WHU-SGG and CNES-GRG 
products respectively, the one cycle inconsistency arises whenever the superposition of 
receiver FCB and satellite FCB exceeds the boundary. The additional iterations of LSM 
and computation of a large matrix demand a long time to finish. Therefore, a fast and 
efficient estimation method is desirable. 


In this contribution, a Kalman filter is employed to speed up the computation. The large 
number of observations are handled epoch by epoch, which significantly reduces the 
dimension of the involved matrix and accelerates the computation. The recursive com- 
putation of the Kalman filter allows real-time applications. As for the one cycle incon- 
sistency, a pre-elimination method is developed to avoid the posterior adjustments and 
iterations of the whole process. The following section describes the theoretical back- 
ground of FCB estimations and Kalman filter. Then a pre-elimination method of one 
cycle inconsistency is proposed following the analysis of temporal stabilities of satellite 
FCBs. A set of FCB products is determined and evaluated by comparing with those of 
CNES and WHU. With the estimated FCBs, the improvements from PPP AR are as- 
sessed. Finally, the methodology is summarized, and an outlook for future research is 
presented. 


2 Methodology 


We start with the basic PPP float mode followed by a description of FCB estimation strat- 
egy. The pre-elimination method of one cycle inconsistency is elaborated as well as an 
introduction to Kalman filter. 


39 


Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution 


2.1 PPP float mode 


In GNSS dual-frequency PPP, the ionospheric-free (IF) combination is routinely em- 
ployed to eliminate the effect of the first-order ionospheric delay. For a satellite s ob- 
served by receiver r, the corresponding pseudorange and carrier phase observation 
equation can be expressed as 


du = pr + dt, — dt* + dT + D, jp — Dir + €p,, (1) 


PF ip = p? + dt, == dt’ + dT + Am (NE ir + Br ir = Bir) + Edif 


where p; indicates the geometric distance between satellite and receiver; dt, and dt? are 
the clock errors of receiver and satellite; dT is the slant tropospheric delay; D, r and Dj; 
are the receiver and satellite specific code hardware delays; Ar and N7;r are the wave- 
length in meters and ambiguity in cycles; B, zr and By, are the receiver-dependent and 
satellite-dependent uncalibrated phase delays; £p,, , €», are the pseudorange and carrier 
phase measurement noise. 


Conventionally, precise orbit and clock products from the IGS analysis center are used 
to remove satellite orbit and clock errors. During the generation of IGS precise products, 
the pseudorange ionospheric-free hardware delay bias Dj; is assimilated into the clock 
offset dt? in accordance with the IGS analysis convention. Due to the fact that pseudor- 
ange measurements provide the reference to clock parameters, the actual receiver clock 
estimate would absorb the ionospheric-free combination of the receiver pseudorange 
hardware delay D, zp. After applying the GNSS precise satellite clock products, equation 
(1) can be rewritten as 


PSip = pÈ dT + dt, Ep, 


= em (2) 
PF ir = pè + dT + dt, + Arr Ns rg + Erp 
where the receiver clock and ambiguity can be re-parameterized as 
dt, = (dt, + Dy tr) (3) 
al E S 
Near = Nr ır + brir — bie (4) 
Dy F 
brir = Brip — - (5) 
IF 
Dir 
bîp = Bir — I (6) 
IF 


—S 
For ambiguity-float PPP solutions, the ionospheric-free ambiguity parameter N, rp is es- 
timated as a real-value constant. Since the estimated ambiguity parameter is a combina- 


40 


Chapter 2 Estimating satellite phase fractional cycle biases based on Kalman filter 


tion of the integer ambiguity, the corresponding code hardware delays, and the uncali- 
brated carrier phase delays at both receiver and satellite ends, the integer property is 
lost. 


2.2 FCB estimation 


For an ambiguity-fixed PPP solution, Noe is usually decomposed into the following 
combination of integer WL and float NL ambiguities for ambiguity fixing 


ch S C 775 
i _ fF Ne wL + fi + fh Ny NL 


—s 
r,IF — 
Air 


(7) 


Note that V7,; is an integer WL ambiguity, which implies that WL FCB will not directly 
contribute to PPP. The WL ambiguity can be resolved by the Hatch-Melbourne-Wüb- 
bena (HMW) combination observable (Hatch 1982; Melbourne 1985b; Wiibbena 1985a), 


fidia hiz  fPrat foPr2 


ae: fit f 


= w (New. + dew — divi) (8) 


where 


faDra Tt f2Dr2 


d = B.., — B, — —— (9) 
oe "t = Awı(fı + f2) 
hDi+ RD 
d$,, = Bj — Bs — ——— (10) 
WE 1 wil + fa) 
=S 
N; wL = New T d, wi, T divı (11) 


If the WL ambiguity is correctly resolved to an integer value based on (8), the NL ambi- 
guity observable can be derived, 


- 


NyNL = N$a T d, NL = dy (12) 
where 
hth 
d, NL = » (Arg By. ip = Dy tr) (13) 
fit+ f 
dy; = = Aue Bir — Dix) (14) 


Equations (11) and (12) serve as the basic model for estimating FCBs. Since they have the 
same structure, a general expression can be formulated as 
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RS =N,—NS = d, — d5 (15) 


for WL and NL linear combinations, respectively. R7 represents the FCB measurements; 
N, denotes the float undifferenced ambiguities; N? denotes the integer part of N,, which 
is the sum of the original integer ambiguity and the integer part of the combined code 
hardware delays and uncalibrated phase delays from both receiver r and satellite s; 
d,and d? denote the receiver and satellite FCBs. 


A set of equations in the form of (15) can be integrated based on a network of reference 
stations. Suppose that there are n satellites tracked at m reference stations, the system of 
equations can be expressed as 


s,1 No = N32 
RA r,1 r,1 1:5 20: eus O 
: : : : H : " : dyi 
= —s,n . . . . Li » : 
Ret Nya -NT 1 = 0 0 ws -1 d 
lee ll (16) 
RS! —s,1 si 0 >œ 1 -1 -- 0 d51 
nm Nym T Nim : N 3 . ; ] : 
: : ; : $ : * : po 
rem | |osn De. OR ee 


sn 
Nm m Nem 


The obtained system of equations is singular. The integer ambiguities N? need to be de- 
termined for each equation on the left side, while one arbitrarily selected FCB should be 
set to zero on the right side. Assuming the float ambiguities N} are precisely estimated, 
N} can be determined by rounding Nj. In this way, the system of equations can be re- 
solved. However, the rounding approach may introduce one cycle inconsistencies. 


2.3 One cycle inconsistency 


Since the FCB estimates are limited in the range of one cycle, e.g. [-0.5, 0.5] and [0, 1] for 
WHU-SGG and CNES-GRGS products respectively, the one cycle inconsistency arises 
whenever the superposition of receiver FCB and satellite FCB exceeds the boundary. A 
simple example is presented to depict the situation. Assuming there is one satellite (d? = 
0.2) tracked at two stations (d, = 0.6 and d, 2 = 0.8), the superposition of FCBs should 
be (R74 = —0.4 and R} = —0.6). However, the R7; would be 0.4 due to the rounding 
approach. As a consequence, the corresponding satellite FCB derived from the two sta- 
tions differs one cycle. The one cycle inconsistency can be detected by examining the 
posterior residuals, as employed in former works. The posterior adjustments are ineffi- 
cient as iterations of the whole process are required. In contrast, we propose to eliminate 
the inconsistency in advance, which shows a clear advantage over posterior adjustment. 
The proposed method consists of the following steps: 
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(1) For all the FCB measurements at a single station, the satellite FCBs are eliminated 
using previous estimates. The underlying assumption is that satellite FCBs are stable 
over successive epochs and can be eliminated to a large extent by previous estimates, 
which will be proved in the next section. 


(2) The residual of each FCB measurement after the subtraction of the satellite FCB yields 
an initial estimate of the receiver FCB. A set of initial receiver FCB estimates could be 
obtained from all simultaneously observed satellites. 


(3) In theory, the receiver FCBs obtained from all simultaneously observed satellites are 
expected to be consistent. Therefore, a one cycle inconsistency can be detected by exam- 
ining the group of initial receiver FCBs. If an inconsistency exists, the corresponding re- 
ceiver FCB will differ by +1 cycle. 


(4) To compensate the one cycle inconsistency in receiver FCBs, the corresponding inte- 
ger ambiguity N‘ is adjusted by +1 cycle. 


(5) The procedures described above can be performed in an iterative way until a con- 
sistent receiver FCB is obtained from all satellite measurements. 


In addition, stations can also be handled individually to obtain clean FCB measurements. 
After the pre-elimination of the one cycle inconsistency, the clean FCB measurements 
are fed to the Kalman filter. 


2.4 Kalman filter 


The Kalman filter addresses the general problem of state estimates of discrete time con- 
trolled processes that are governed by a linear stochastic difference equation (Kalman 
1960). The theory has been well studied and widely applied (Yang 2006; Yang et al. 2010). 
Since the Kalman filter is based upon the theory of least squares, it is theoretically pos- 
sible to calculate the same solution as LSM and versa vice. However, the design and 
normal equation matrix will be huge in LSM, considering the large number of observa- 
tions from hundreds of reference stations and thousands of epochs. Computation of the 
large matrix is time-consuming, which is inefficient and unfavorable for real-time appli- 
cations. In contrast, the large number of observations are handled epoch-by-epoch in a 
Kalman filter, which significantly reduces the dimension of the involved matrix and ac- 
celerates the computation. One additional advantage of Kalman filter is its real-time ca- 
pability. This is of particular interest for real-time PPP AR and its applications. There- 
fore, a Kalman filter is employed in this work. The design and flowchart of a Kalman 
filter are depicted in Figure 1. In our case, Rẹ is the input FCB measurement for the Kal- 
man filter, while d, and d? are the output estimates. 
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Discrete linear system 
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Figure 1 Flow chart of the proposed Kalman filter. The dashed blocks represent additional steps adopted 
in this study 


Two additional steps are added to the standard Kalman filter. The first step aims at es- 
tablishing the dynamic model and determining the state transition matrix by analyzing 
the temporal stability of existing FCB products. Satellite WL FCBs are stable over several 
days, which can be characterized as a constant parameter on a daily basis. Satellite NL 
FCBs are considered as constant within 15 minutes but exhibit small variations over 15- 
minute intervals. Therefore, NL FCBs are characterized as random walk processes. The 
second step is introduced to eliminate the one cycle inconsistency of FCB measurements. 
Only clean measurements are sent to the Kalman filter to avoid iterations. Note that a 
constraint is imposed on the Kalman filter to eliminate the rank deficiency. A satellite 
FCB is selected and set to zero, which is accomplished by adding a pseudo observational 
error equation with zero variance (Yang et al. 2010). 


3 Results, comparison, and analysis 


In this section, the temporal stability of FCB is first analyzed based on existing products. 
Then, the proposed pre-elimination method of one cycle inconsistency is demonstrated. 
Ten sets of FCBs are computed and evaluated with respect to existing products. At last, 
with the estimated FCBs, the improvements from PPP AR are assessed. 
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3.1 Temporal stability of FCB 


In order to study the characteristics of satellite FCBs, products from CNES GRG and 
WHU SGG are employed. Both organizations provide WL FCBs derived from the same 
strategy, while NL FCBs are only available for SGG products since they are absorbed by 
satellite phase clock offsets in GRG products. The products for entire 2016 are down- 
loaded and analyzed in this study. Due to the rank deficiency of FCB equations, a re- 
ceiver FCB for GRG products and a satellite FCB for SGG products are selected as refer- 
ences. Therefore, a single differencing process across satellites is performed to remove 
the datum before comparison. The single-differenced WL FCBs are presented in Figure 
2. Note that the two products have opposite signs. 


SD-WL-FCB [Cycle] 


SD-WL-FCB [Cycle] 


50 100 150 200 250 300 350 
Day of Year 


Figure 2 Single-differenced WL FCBs of CNES GRG (top) and WHU SGG (bottom) products in 2016. The 
lines represent individual satellite WL FCBs with respect to that of G01 


The WL FCBs of satellites are remarkably stable over days with standard deviations 
(STDs) of less than 0.08 cycles for both products. An exception is observed for G32 satel- 
lite. The STD of G32 is significantly larger than the others with 0.2 cycles. The reason is 
that a new Block IIF satellite (G32) was launched on 5 February 2016 and activated since 
10 February. After the replacement, the WL FCB exhibits a similar temporal stability as 
the other satellites. The statistics further confirm the conclusion that the satellite WL FCB 
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is stable over at least several days. Therefore, WL FCBs are estimated as constant on a 
daily basis in the proposed Kalman filter. 


NL FCB products are estimated with respect to specific IGS precise clock products in the 
WL/NL FCB based method. Due to the daily boundaries of precise clock products, NL 
FCB estimates are only continuous within one day. The single-differenced SGG NL FCBs 
on DoY 001, 2016, are presented in Figure 3. The top panel shows the raw SD NL FCBs 
for each satellite while G01 is selected as reference. It can be seen that the datum changes 
frequently, which introduces a difference of one cycle between two consecutive epochs. 
After adjustments, the NL FCB series are continuous, as presented in the middle panel. 


One can easily discern that satellite NL FCBs are also stable. The STDs of all single-dif- 
ferenced satellite NL FCBs are within + 0.05 cycles for DoY 001. In addition, daily STDs 
of all satellite NL FCBs for the whole year 2016 are calculated. The bottom panel depicts 
the distribution of daily STDs. 98.37% of the daily STDs are below 0.15 cycles while 
97.02% are below 0.1 cycles. These results show that satellite NL FCBs may be more sta- 
ble than reported in former research. Note that a small number of daily STDs may reach 
0.5 cycles. These abnormal STDs may be caused by maneuvers during satellite eclipse or 
unmodeled errors. And we could not exclude the possibility that there may be blunders. 
Nevertheless, satellite NL FCBs are modeled as random walk processes in the proposed 
Kalman filter. Also note that the characteristics of receiver NL FCB depend on the re- 
ceiver type, environment, and other factors. It is hard to characterize receiver NL FCBs 
with a general model. Therefore, they are modeled as white noise. The temporal stability 
of satellite FCBs also implies that it can be utilized for pre-elimination of the one cycle 
inconsistency. 
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Figure 3 Raw (top) and adjusted (middle) Single-differenced NL FCBs of SGG products on DoY 001. The 
bottom panel presents the histogram of all daily STDs in 2016. 


3.2 Pre-elimination of the one cycle inconsistency 


For a successful pre-elimination of the one cycle inconsistency, the satellite FCBs should 
be removed in advance. The key issue is if or at what extent the satellite FCBs of the 
current epoch can be counteracted by the previous estimation. In the proposed pre-elim- 
ination method, the differences of satellite FCBs between two successive epochs will be 
absorbed by receiver FCB. Therefore, epoch differences of satellite FCBs should be small 
enough to avoid impacts on the detection of the one cycle inconsistency. In order to val- 
idate the hypothesis, an additional epoch differencing process is performed on the sin- 
gle-differenced FCBs. The histograms of epoch differences for WL FCBs in 2016 are pre- 
sented in Figure 4. 
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Figure 4 Histograms of epoch differenced WL FCBs of GRG (top) and SGG (bottom) products in 2016 


99.78% and 95.8176 of the differences are within + 0.05 cycles for GRG and SGG products 
respectively. The max differences, with a magnitude of 0.386 cycles, are found in SGG 
products. Since these values are much smaller than one cycle, it will scarcely affect the 
detection and elimination of the one cycle inconsistency. Note that the distribution of 
epoch differences for GRG products is more concentrated around zero than that of SGG, 
which may indicate better quality. 


The epoch differencing process is also conducted for NL FCB products, as shown in Fig- 
ure 5. 99.73% of the differences are within + 0.05 cycles for SGG products. The result is 
even better than that of WL FCB. This can be attributed to the fact that the intervals (15 
minutes) of NL FCB are much shorter than those of WL FCB (24 hours). 
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Figure 5 Histogram of epoch differenced NL FCBs of SGG products 
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Figure 6 Pre-elimination of the one cycle inconsistency for station IOOE on 2016/02/26. The top panels rep- 
resent receiver WL FCB adjustment, while the bottom ones represent that of receiver NL FCB. The left and 
right panels represent the raw and adjusted measurements respectively 


According to the above analysis of satellite FCB products, it can be concluded that the 
major parts of both WL and NL satellite FCBs can be mitigated by previous estimations. 
After the removal of satellite FCBs, the one cycle inconsistency among FCB measure- 
ments can be detected by examining the residuals. After the detection of inconsistency, 
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it is eliminated by adjusting the corresponding integer ambiguity. An example of the 
process is shown in Figure 6. One can easily observe that the one cycle inconsistencies 
occur on several satellite measurements and can be effectively eliminated by the pro- 
posed method. It should be noted that elimination of the one cycle inconsistency will not 
improve the precision of FCB estimation as the fractional part remains the same. In this 
sense, the magnitude of the epoch difference is less of concern as long as it is sufficient 
to detect potential inconsistencies. Also, note that this can be accomplished by posterior 
residuals adjustments. However, pre-elimination shows a clear advantage over posterior 
adjustment in terms of efficiency. 


In this way, all epochs and all station data can be processed individually to obtain clean 
FCB measurements. For post-processing, the above procedure can be simplified. Daily 
WL FCB measurements can be adjusted simultaneously as the receiver WL FCBs are also 
constant within one day. For the initialization of the Kalman filter, a preliminary set of 
satellite FCBs should be provided for the first epoch. The initial satellite WL FCBs are 
taken from the estimations of the last day, while the initial satellite NL FCBs are deter- 
mined with the route method (Li et al. 2015a). 


3.3 Comparison of FCB products 


A globally distributed network consisting of about 200 stations is selected from IGS net- 
work, as shown in Figure 7. Observations from Day of Year (DoY) 52 to 61 (2016/02/21 - 
2016/03/01) are processed to determine GPS satellite FCBs. GFZ final precise products 
are applied to remove satellite orbit and clock errors. The other errors are corrected ac- 
cording to the IGS standard error models (Kouba 2009). FCB estimations are conducted 
in three sequential steps. First, float ambiguities containing FCBs are obtained by HMW 
combinations and PPP processing from the network of reference stations. Secondly, FCB 
measurements are generated from the float ambiguities. Thirdly, the proposed Kalman 
filter is adopted to estimate FCBs from all FCB measurements. 


The quality of FCB estimation can be indicated by the posterior residuals of FCB meas- 
urements. In general, a highly consistent FCB estimation is expected if the residuals are 
close to zero. Figure 8 shows the distribution of the residuals of FCB estimations for the 
ten days. It can be seen that both histograms are symmetric and nearly centered at zero. 
The Root Mean Square (RMS) of WL residuals is 0.09 cycles while the RMS of NL resid- 
uals is about 0.05 cycles, which indicates a good consistency between the estimated FCBs 
and the input float ambiguities. The total number of input float ambiguities is 934908. 
Figure 9 shows the usage rates of WL and NL float ambiguities for each satellite. The 
averaged usage rate for WL is 96.66%, while that for NL is 98.17%. We find that the qual- 
ity of NL FCBs is better than that of WL. The possible reason is that, WL FCBs are esti- 
mated as daily invariants and easily affected by pseudorange noise, while the NL FCBs 
are updated every 15 minutes and derived from much more precise carrier phase meas- 
urements. 
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For better visualization, the derived satellite FCBs are shifted with integer cycles, as pre- 
sented in Figure 10. It can be seen that the derived satellite FCBs show similar temporal 
stability as in the above analysis. Satellite WL FCBs exhibit extremely small variations 
over the ten days period, while NL FCB shows a larger variation over tens of minutes. 


Figure 7 Geographical distribution of the selected reference stations 
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Figure 8 Distributions of posterior residuals after WL (top) and NL (bottom) FCB estimation 
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Figure 9 Usage rates of the float WL and NL ambiguities 
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Figure 10 Time series of satellite WL (top) and NL (bottom) FCBs for all the other 31 satellites while G01 is 
selected as reference 
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In order to validate the external accuracy of our estimation, the derived satellite FCBs 
are also compared with those of SGG and GRG. As discussed before, WL FCB measure- 
ments are obtained from HMW combinations and are relatively independent of PPP pro- 
cessing. The determined WL FCBs should be consistent across all products. Figure 11 
presents the comparison of the derived satellite WL FCBs with those of GRG and SGG 
products. 
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Figure 11 Histograms of the differences between our WL FCB estimation and SGG (top) and GRG prod- 
ucts (bottom) 


It can be seen that our WL FCB estimates show good consistency with those of SGG and 
GRG. All of the differences are within + 0.1 cycles. Compared with SGG products, 76.7% 
of the differences are within + 0.05 cycles with an RMS of 0.05 cycles. Compared with 
GRG products, 91.0% of the differences are within + 0.05 cycles with an RMS of differ- 
ences of 0.03 cycles. Based on the above analysis, we can safely conclude that there is no 
systematic bias between our WL FCB estimates and those of SGG and GRG. The differ- 
ences are actually minimal and negligible. However, one can realize that our WL FCB 
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estimates agree better with GRG products than those of SGG. Note that the RMS of dif- 
ferences between SGG and GRG is 0.03 cycles. We suspect that the GRG products may 
perform slightly better in practice. 
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Figure 12 Histogram of the differences between our NL FCB estimation and SGG products 


NL FCBs are only available for SGG products. The derived satellite NL FCBs are com- 
pared with SGG products, as shown in Figure 12. It can be seen that our NL FCBs agree 
well with SGG NL FCB products. 97.9 % of the differences are within + 0.1 cycles while 
70.4% of the differences are within + 0.05 cycles. The RMS of the differences is 0.05 cycles, 
corresponding to 5.1 mm. The discrepancy between the two results may be ascribed to 
different PPP processing strategies. Since NL FCB measurements are directly derived 
from PPP float ambiguities, discrepancies between error correction models, such as trop- 
ospheric models, may introduce small systematic biases. Another possible reason could 
be the different distributions of reference stations. Since FCB estimates are contaminated 
by unknown temporally and spatially correlated errors, these unknown errors may 
change with the distribution of reference stations. GRG products should be employed 
for comparison in further investigations. Note that the discrepancy will not affect PPP 
AR at user end if the same error model as at the server end is used. 


3.4 PPP AR solution 


In order to validate our FCB estimates and to assess the improvement by ambiguity res- 
olution, all IGS network stations are processed in PPP float and AR mode with the ob- 
tained satellite FCB estimates. The 24 hours observations from over 500 stations in the 
ten days are divided into three-hour sessions. The solutions with data integrity less than 
80 percent, wrongly fixed, and incomplete convergence are removed. After the pre-pro- 
cessing, there are 22953 sessions. The performance in terms of convergence time and 
positional accuracy is evaluated under different confidence levels for the sake of relia- 
bility (Lou et al. 2016). 
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Figure 13 Convergence performance of GPS PPP float and AR solutions based on 22953 three-hour ses- 
sions under 95% confidence levels 


Table 1 Accuracy comparison of GPS PPP float and AR solutions based on 22953 three-hour sessions un- 
der different confidence level [Unit: cm] 


Time Solut. 
50% 68% 95% 50% 68% 95% 50% 68% 95% 


Float 157 251 729 0.71 1.09 2.67 2.06 3.16 8.01 
One-hour AR 0.43 0.67 2.68 0.37 058 1.62 145 226 5.87 
Improv. 73% 73% 63% 48% 47% 39% 30% 28% 27% 

Float 0.95 149 397 044 067 157 141 211 48 

Two-hour AR 0.37 0.57 148 029 045 112 1.09 167 39 
Improv. 61% 62% 63% 34% 33% 29% 23% 21% 19% 
Float 073 11 253 034 051 116 121 1.6 3.76 

Three-hour AR 0.34 052 133 026 039 0.96 0.91 1.39 3.29 
Improv. 53% 53% 47% 24% 24% 17% 25% 21% 13% 


Figure 13 presents the positional errors of GPS PPP float and AR solutions based on the 
statistics over all the sessions. It can be seen that the convergence time is significantly 
shortened by ambiguity resolutions, especially for the east and up components. It takes 
64 minutes for float solutions to converge to three-dimensional 10 cm accuracy, while 
that for AR solutions is only 48.5 minutes, corresponding to an improvement of 24%. In 
order to assess the improvements with respect to time length, the positional errors are 
calculated for one-hour, two-hour, and three-hour solutions, as presented in Table 1. It 
can be seen that PPP ambiguity resolutions are able to enhance the accuracy for all the 
schemes. The most significant improvement is found for the east component. Due to the 
design of satellite navigation systems, the accuracy of the east component is worse than 
that of north component in mid-latitude because of weaker model strength. Ambiguity 
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resolutions can improve this situation by imposing tight constraints on ambiguity pa- 
rameters. However, the improvements decrease as the time length increases. The im- 
provements of ENU components for one-hour solutions are (70%, 45%, 28%), while that 
for two-hour solutions and three-hour solutions are (62%, 32%, 21%) and (51%, 22%, 
20%), respectively. The reason is that the model strength of the float solution is improved 
with more observations. Also note that the improvements differ with respect to the con- 
fidence level. They are less significant under 95% level, which implies the PPP AR is not 
effective for some stations. Since FCB estimates are contaminated by temporally and spa- 
tially correlated errors, a small and dense network is preferred for better performance. 


4 Conclusions 


In this contribution, estimating satellite FCBs based on the Kalman filter is proposed and 
demonstrated. Since the Kalman filter is based upon least squares methods (LSM), it is 
theoretically possible to calculate the same solution as for the commonly used LSM. In 
the proposed Kalman filter, the large number of observations is handled epoch-by- 
epoch, which significantly reduces the dimension of the involved matrix and accelerates 
the computation. Hence it outperforms the commonly used LSM in terms of efficiency. 
In order to avoid iterations caused by the one cycle inconsistency among FCB measure- 
ments, a pre-elimination method is developed based on the temporal stability of satellite 
FCBs. The pre-elimination method shows a clear advantage over post-residual adjust- 
ment, which further improves the efficiency. 


A globally distributed network consisting of about 200 IGS stations has been selected to 
determine GPS satellite FCBs. The estimated WL FCBs have a good consistency with 
existing WL FCB products (e.g., CNES-GRG, WHU-SGG). The RMS of differences with 
respect to GRG and SGG products are 0.03 and 0.05 cycles, which indicates the con- 
sistency of the proposed approach. For satellite NL FCB estimates, 97.9% of the differ- 
ences with respect to SGG products are within + 0.1 cycles. The RMS of the differences 
is 0.05 cycles. These results prove the efficiency of the proposed approach. 


The state-based approach of the Kalman filter also allows for more realistic modeling of 
stochastic parameters, which will be investigated in future research. 
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Abstract: Due to the rapid deployment of the Galileo constellation, Galileo is now able to contrib- 
ute to GNSS precise point positioning (PPP) ambiguity resolution (AR) with 17 operational satel- 
lites as of December 2017. We estimate the satellite fractional cycle bias (FCB) based on globally 
distributed MGEX stations and assess the Galileo FCB quality by a comparison with that of GPS 
and BDS. Results of 60 days indicate that the quality of Galileo wide-lane (WL) FCB is better than 
GPS and BDS in terms of data usage rate, residual distribution, as well as standard deviation of 
daily estimates. The RMS of Galileo WL FCB residuals is 0.071 cycles, while that of GPS and BDS 
are 0.089 and 0.117 cycles respectively. The standard deviation of Galileo daily WL FCB is 0.010 
cycles, while that of GPS and BDS are 0.018 and 0.043 cycles. We attribute the better quality of 
Galileo WL FCB to its signal modulation, AltBOC, which significantly compresses the multipath 
effect for pseudorange measurement. Within the Galileo constellation, the performance of In-Or- 
bit Validation (IOV) satellites WL FCB is worse than that of Full Operational Capability (FOC) 
satellites as a result of a reduction in the power of the transmitted signal. The performance of the 
two highly eccentric satellites is comparable to other FOC satellites. The overall quality of Galileo 
narrow-lane (NL) FCB is slightly worse than that of GPS but better than that of BDS. The RMS of 
Galileo NL FCB residuals is 0.062 cycles, while that for GPS and BDS are 0.050 and 0.086 cycles 
respectively. In addition, the NL FCB quality of FOC, IOV (except E19), as well as the two eccen- 
tric satellites, shows no significant difference in terms of data usage rates and residuals. Galileo 
PPP AR solutions are conducted at 20 MGEX stations with three-hour sessions for ten days. The 
positional biases of AR solutions are 0.7, 0.6 and 2.1 cm for east, north and up components re- 
spectively, while those for float solutions are 2.1, 1.1 and 2.7 cm, corresponding to improvements 
of 67, 45 and 22% respectively. These results demonstrate that currently Galileo FCB can be esti- 
mated with an accuracy comparable with GPS and BDS, and the Galileo observations can bring 
an obvious benefit to ambiguity-fixed PPP. 
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1 Introduction 


The fact that double differenced ambiguities in network solutions can be resolved to in- 
teger values lays the foundation of integer ambiguity resolution for precise point posi- 
tioning (PPP). By estimating the fractional cycle biases (FCB) at the server and applying 
them to single differenced PPP at the user, the PPP integer ambiguity resolution can be 
achieved, convergence time significantly shortened, and the accuracy improved (Teunis- 
sen and Khodabandeh 2015). Ge et al. (2008) proposed to estimate the FCBs of the single 
differenced ambiguities in wide-lane (WL) and narrow-lane (NL) form by averaging the 
fractional parts of all involved float single-differenced WL and NL ambiguity estimates. 
Instead of an averaging process, a least squares method in an integrated adjustment was 
adopted to enhance the FCB estimates in Li et al. (2015a) and Li and Zhang (2012). Dif- 
ferent from the above WL/NL FCB approaches, Collins et al. (2008) developed a decou- 
pled clock model by separating satellite clocks for code and phase observations. Simi- 
larly, an integer phase clock model in which the NL FCBs were assimilated into receiver 
and satellite clock estimates of a network solution was developed (Laurichesse et al. 
2009; Loyer et al. 2012). These PPP ambiguity resolution (AR) techniques have been 
proven equivalent in theory (Shi and Gao 2014), and the positional biases have been 
demonstrated to be minimal (Geng et al. 2010). 


In addition to the development of functional models, PPP AR has been extended to mul- 
tiple GNSSs. For GLONASS, special care has to be taken to deal with the receiver inter- 
frequency bias (IFB) between satellites originating from the frequency division multiple 
access strategies. Usually, it is achieved by using a network of homogeneous receivers 
(Geng and Shi 2016; Liu et al. 2017b; Yi et al. 2017). Also, BDS-2 code measurements 
suffer from the satellite induced pseudorange variations which are elevation- and fre- 
quency-dependent. Correction models based on multipath combination have been 
proven effective (Wanninger and Beer 2015). Currently, only regional or IGSO/MEO PPP 
AR is achievable for BDS (Gu et al. 2015a). Tegedor et al. (2016) initially assessed the 
Galileo PPP AR with four In-Orbit Validation (IOV) satellites using regional stations 
around Europe. Li et al. (2017b) estimate the FCBs for the four global systems and assess 
the performance of combined PPP AR. 


Compared with GPS and BDS, the study of Galileo PPP AR has been limited by the num- 
ber of available satellites. Galileo is currently at the initial service phase. After the de- 
commission of Galileo In-Orbit Validation Element-A (GIOVE-A) and GIOVE-B, four 
IOV spacecraft were launched between 2011 and 2012. Two Full Operational Capability 
(FOC) Galileo satellites were launched, however, into wrong, highly eccentric orbits in 
2014 (Sonica et al. 2017). While the two satellites are unlikely to ever become a part of 
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the operational constellation, they offer proper navigation signals and broadcast naviga- 
tion messages, which allows for real-time navigation and PPP (Montenbruck et al. 2017). 
In 2015 another six, in May 2016 another two, and in Nov 2016 another four FOC Galileo 
satellites were successfully launched. The full constellation with 30 satellites is expected 
to be realized by 2020. As of December 2017, there are 17 active Galileo satellites as 
shown in Table 1. Note that E20 is excluded as it transmits only E1 signal and is not 
available for dual-frequency positioning. These satellites are separated into three groups: 
four IOV, two highly eccentric (ECC) and eleven FOC satellites. Galileo provides an open 
signal in the E1 band and a wideband signal covering the E5 a&b band. The Alternating 
Binary Offset Carrier (AltBOC) modulation can either be tracked as a composite signal 
or as distinct signals in the E5a and E5b sub-bands. Zaminpardaz and Teunissen (2017) 
compared the signal power, multipath performance, code and phase noise between IOV 
and FOC satellites. Results show that their characteristics differ significantly. It is neces- 
sary to investigate the impact on FCB estimation. In addition, it was found that the orbit 
quality of the two highly eccentric satellites is comparable with that of IOV and FOC 
satellites (Sonica et al. 2017; Steigenberger and Montenbruck 2017). Thus, it is also worth- 
while to evaluate the performance of FCB estimation of the two eccentric satellites. 


Table 1 List of Galileo satellites as of December 2017 (https://www.gsc-europa.eu/system-status/Constella- 
tion-Information) 


Satel- PRN Status 
GIOVE A,B Retired 
E11, E12, E19 Operational 
IOV 
E20 Unavailable 
E18, E14 Launched into incorrect orbit, denoted as “ECC” 
satellites 
FOC E01, E02, E03, E04, E05, E07, E08, E09, E22, . 
Operational 
E24, E26, E30 
E21, E25, E27, E31 Commissioning 


In the context of multi-GNSS PPP, the aim of this study is to extend the PPP ambiguity 
resolution to Galileo based on the globally distributed MGEX observations. Addition- 
ally, how FCB estimation can benefit from the advanced technique of Galileo signal mod- 
ulation is investigated. Furthermore, within the Galileo satellite constellation, we also 
conduct a comparison among IOV, highly eccentric and FOC satellites. With the esti- 
mated FCBs, the improvements from Galileo PPP AR are evaluated. Finally, the meth- 
odology is summarized, and an outlook for future research is presented. 
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2 Methodology 


We start with a presentation of the basic PPP float mode, followed by a description of 
PPP data processing for ambiguity resolution. Then our FCB estimation strategy is in- 
troduced. With the obtained FCB products, the PPP AR algorithm at the user end is elab- 
orated. 


2.1 Ambiguity-float PPP mode 


In GNSS dual-frequency PPP, the ionospheric-free (IF) combination is routinely em- 
ployed to eliminate the effect of the first-order ionospheric delay. For a satellite s ob- 
served by receiver r, the corresponding pseudorange and carrier phase observation 
equations can be expressed as 


Pa = p; + dt, — dt” + dT + Dy, jp — Dig + €p,, (1) 


PF ip = p + dt, xd dts + dT + Am (NE ir + Br ir = Bir) + Edif 


where p; indicates the geometric distance between satellite and receiver; dt, and dt? are 
the clock errors of receiver and satellite; dT is the slant tropospheric delay; D, r and Dj; 
are the receiver and satellite specific code hardware delays; Ar and N7; are the wave- 
length in meter and integer ambiguity in cycle; B, ¡r and Bfp are the receiver-dependent 
and satellite-dependent uncalibrated phase delays; £p,, , £»,, are the pseudorange and 
carrier phase measurement noise. 

Conventionally, precise orbit and clock products from the IGS analysis center are used 
to remove satellite orbit and clock errors. The pseudorange ionospheric-free hardware 
delay bias Dj; is assimilated into the clock offset dt? in accordance with the IGS analysis 
convention. Due to the fact that pseudorange measurements provide the reference to 
clock parameters, the actual receiver clock estimate would absorb the ionospheric-free 
combination of the receiver pseudorange hardware delay D, r. After applying the GNSS 
precise satellite clock products, equation (1) can be rewritten as 


Pj = ps + dT + dt, + ep, 


E = (2) 
Dir = p; + aT + dt, + DEN T t Edif 
where the receiver clock and ambiguity can be re-parameterized as 
dt, = (dt, + Dy tr) (3) 
7 E S 

Nyair = Near + bye — Dir (4) 

Dr iF 
brir = Brr — - (5) 

IF 


62 


Chapter 3 Estimating and assessing Galileo satellite fractional cycle bias for PPP ambiguity resolution 


S S Dij F 
bjr = Bir = p (6) 
IF 
Since the estimated ambiguity parameter is a combination of the integer ambiguity, the 
corresponding code hardware delays, and the uncalibrated carrier phase delays at both 
receiver and satellite ends, the integer property is lost. 


2.2 FCB estimation 


In order to fix PPP ambiguity, the term N, ir is usually decomposed into the following 
combination 


ch 


S C TS 
Ls f = f Nr wL + fi + h N; NL 
N; iF = GE | os (7) 
IF 

where Nw, is the integer WL ambiguity and N, nz is the float NL ambiguity. Usually, 
the WL ambiguity is resolved by the Hatch-Melbourne-Wiibbena (HMW) combination 
observable (Hatch 1982; Melbourne 1985b; Wübbena 1985a). With the fixed WL ambigu- 
ity, the float NL ambiguity can be derived and tested whether it is also fixable. An iono- 
spheric-free ambiguity is fixed only when both its WL and NL ambiguities are fixed. 


The float WL ambiguity can be derived from 


fii. fi, fiPfa + Ps, 


—s 
HMW = —ı = 28,577 N 
h-fe h tfe TOUS 


(8) 


where Ay wı = rare is the wavelength of WL ambiguity. The float WL ambiguity can be 
1 72 


further decomposed as 


ELS 


Ny wL = New + d, wL Ex diy, (9) 


fiDr1 + f2Dr.2 
Awi + fa) 


fi Dj + f2D3 
wL + f2) 


dy,wı E Bra = By? = (10) 


diy, = Bi — B? — (11) 
HMW combinations eliminate the geometric distance, satellite and receiver clock errors, 
as well as atmospheric delays. The key factors that affect HMW accuracy would be code 
noise and multipath effect. The noises of code measurements will be smoothed by the 
more precise phase measurements when averaging the HMW for a continuous arc. 
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If the WL ambiguity N/y,, is correctly resolved to an integer value and the float IF am- 


biguity N; ır is obtained from PPP solution, the float NL ambiguity observable N; y; can 
be derived based on (7) 


xs fi th fe 


Ny y, = 2— PAN, rIF  R- he = Newt = Nea + deni — Ani (12) 

where 
d, nL = mer the —— (Arg B;, IF — D,.1r) (13) 
dy, = S the (AırBir — Dig) (14) 


As we can see, there are no code measurements directly involved in the calculation of 
NL ambiguities. Since the phase measurement is very accurate at millimeter level, the 
unmodeled errors in PPP will be important for the quality of NL ambiguity. 

Equations (9) and (12) serve as the basic model for estimating FCBs. Since they have the 
same structure, a general expression can be formulated as 


RS = N, — Ñ; = d, — d5 (15) 


for WL and NL linear combinations, respectively. R7 represents the FCB measurements; 
N, denotes the float undifferenced ambiguities; V? denotes the integer part of N,, which 
is the sum of the original integer ambiguity and the integer part of the combined code 
hardware delays and uncalibrated phase delays from both receiver r and satellite s; d, 
and d? denote the receiver and satellite FCBs. Considering the temporal stabilities of 
code hardware delays and uncalibrated carrier phase delays, the satellite FCBs can be 
estimated as constant values for specific time intervals, e.g., one day for WL FCB and 15 
minutes for NL FCB (Ge et al. 2008). 


A set of equations in the form of (15) can be generated based on a network of reference 
stations. Suppose that there are n satellites tracked at m reference stations, the system of 
equations can be expressed as 


sl x = N32 
Ry T,1 r,1 1 0 —1 
: n u d; 1 
: —sn a : : : 
Ry Nya — Net 1 = 0 0 =» -1 d 
Az =|) 43 or os 4 FIDES (16) 
R31 —s,1 csi 0 1 -1 d51 
m Nym = Nr m . . : 
: : : : : : u : an 
RES —sn Ow 10 w -1 


sn 
Nym — Nem 
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The obtained system of equations is singular. Assuming the float ambiguities N? are pre- 
cisely estimated, N} can be determined by rounding N; and one arbitrarily selected FCB 
should be set to zero. In this way, the system of equations can be resolved. 


2.3 PPP AR at the user end 


With the above-estimated FCBs, the satellite FCBs can be corrected, and the single-dif- 
ferenced PPP ambiguities are aimed to be fixed. It is noted that, although single differ- 
encing operations could be conducted on raw observation equations and float ambigui- 
ties level, single differencing at the ambiguity level may be better in practice as it 
provides more flexibility of the choice of a reference satellite and the formation of satel- 
lite pairs. Therefore, we recommend performing the single differencing at ambiguity 
level, 


N,-N, = N?" + asm (17) 


For WL float ambiguities, they can be directly fixed by the rounding approach (Ge et al. 
2008). For NL float ambiguities, they are fed into the LAMBDA algorithm (Teunissen et 
al. 1997). If not all the float NL ambiguities can be fixed by the LAMBDA method, partial 
ambiguity resolution can be employed (Li and Zhang 2015; Teunissen et al. 1999). 

For sequential PPP AR processing realized with Kalman filter, such as real-time applica- 
tions, there exists a discrepancy concerning whether the states of the fixed solution 
should be used for the time update of the subsequent epoch. The float and fixed solution 
loops are independent. Even though the fixed solution is obtained, the states of the float- 
ing loop do not change accordingly (Wang et al. 2017). We argue that the benefits of 
ambiguity resolution are not fully exploited in this approach. On the other hand, the 
fixed ambiguity of the current epoch should be tightly constrained to the succeeding 
epoch in a continuous arc (Takasu and Yasuda 2010). This approach utilizes all the states 
of the fixed solution to the next time update of the float solution and is introduced for 
PPP AR here. Once the WL and NL integer ambiguities are resolved and validated, a 
tight constraint can be reconstructed and imposed on the ionospheric-free float ambigu- 
ities 


0=N,-N, —-NS+N"—-d5™, oy (18) 


In this way, a fixed solution can be obtained at the user. It is noted that only the NL FCB 
is used in the reconstructed fixed ambiguity and has direct influence on the estimated 
parameters. In this sense, the NL FCB has to be estimated as precisely as possible. 
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3 Data and processing strategy 


The International GNSS Service (IGS) established the Multi-GNSS Experiment (MGEX) 
in order to prepare operational services for new and upcoming GNSS, such as the Euro- 
pean Galileo, Chinese BeiDou, and regional systems, such as Japanese QZSS and Indian 
NAVIC system (Montenbruck et al. 2017). Currently, there are over 210 MGEX stations. 
Figure 1 shows the geographic distribution of all the MGEX stations. Among these sta- 
tions, 208 stations support Galileo signal tracking, and 165 stations support BDS signal 
tracking. The MGEX stations, equipped with multiple brands of professional receivers, 
provide almost full and continuous tracking of Galileo dual-frequency signals from the 
17 satellites. Therefore, all the data collected by the MGEX network from June 4 to Au- 
gust 2, 2017, in total 60 days, are processed. 


Figure 1 Geographic distribution of all the MGEX stations. Almost all stations support Galileo, while the 
number of BDS stations, represented by the red circles, is slightly smaller. 


In the processing, E1/E5a are used for Galileo, while L1/L2 and B1/B2 are used for GPS 
and BDS respectively. These signals were chosen according to a high availability, and in 
accordance with the principle of orbit and clock generation (Uhlemann et al. 2016). In 
order to obtain accurate float ambiguities, a GPS/Galileo/BDS combined PPP mode is 
used at the server end. Inter-system biases are estimated per station and per system. The 
combination of multiple GNSSs increases the data usage rates of the incomplete Galileo 
and BDS systems, as the number of observable satellites for either of the two systems 
may be not enough for a standalone PPP solution (Li et al. 2017a). Figure 2 shows the 
global distribution of the number of observable Galileo satellites with a cut-off elevation 
of 10 degrees. It can be seen that a standalone Galileo PPP solution is not feasible in a 
few regions. It is noted that this phenomenon is temporary, and significant improve- 
ments can be expected when the constellation is completed. In addition, multi-GNSS 
PPP is robust against outliers and improves the precision of the float ambiguities. A for- 
ward and backward filter is used to avoid the convergences of ambiguity parameters. 
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For the sake of comparison, the basic settings of data edit are the same for all three sys- 
tems: the cut-off elevation angle is set to 10 degrees; float ambiguities with an elevation 
below 30 degrees or with standard deviation (STD) larger than 0.1m are removed. It is 
noted that the BDS satellite-induced code multipath effects are corrected for IGSO and 
MEO satellites according to Wanninger and Beer (2015), while GEO satellites are ex- 
cluded from the processing. The third generation BDS satellite, which are no longer af- 
fected by such effects (Lei et al. 2017), are also excluded since the data is not publicly 
available. 


3 - 5 6 7 8 


Figure 2 Distribution of the number of observable Galileo satellites at DoY 160, 2017 with a cut-off eleva- 
tion of 10 degrees 


Throughout the processing, MGEX precise products, including precise satellite orbit, 
clock, and earth rotation parameters, provided by GFZ (Uhlemann et al. 2016) are used. 
The satellite phase center offsets and variations are corrected according to the IGS an- 
tenna file (igs14_1949.atx). As for the receiver antenna phase center offsets and varia- 
tions, the correction values for GPS are employed for both Galileo and BDS in accordance 
with the principle of orbit and clock generation (Prange et al. 2017). The detail of our 
FCB estimating strategy could be found in Xiao et al. (2018). 


4 Results, comparisons, and discussions 
In this section, the quality of Galileo WL/NL FCB estimation is thoroughly analyzed. We 
further compared the FCB estimations with that of GPS and BDS in order to characterize 


the Galileo signal properties. The generated FCBs are used for static and kinematic PPP 
AR solutions. 
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4.1 Satellite WL FCB estimation 


Figure 3 presents the time series of Galileo satellite WL FCBs during the 60 days. For 
better visualization, the derived satellite FCBs are shifted with integer cycles. It can be 
seen that the derived Galileo satellite WL FCBs are quite stable over time, with a STD of 
0.010 cycles. It would be sufficient to broadcast daily WL FCB, which is similar to other 
GNSSs (Ge et al. 2008). For comparison, we also present the STD statistics for GPS and 
BDS in Figure 4. Note that BDS C07 and C11 are excluded from the statistics as there 
exist jumps in the time series due to unknown reasons, which has been reported by Li et 
al. (2017a) and Li et al. (2017b). The average STDs are calculated for GPS and BDS, while 
an individual STD is calculated for each Galileo satellite. The STDs of WL FCB are 0.018 
cycles for GPS and 0.041 cycles for BDS. This implies that Galileo WL FCB is less noisy 
than GPS and BDS. We attribute the better performance of Galileo WL FCBs to its signal 
modulation. The advanced AItBOC modulation can significantly compress the multi- 
path effect for pseudorange measurements and decrease the fluctuations of HMW com- 
binations, especially for low elevations (see Figure 5). Also, we notice that there is a dif- 
ference between FOC and IOV satellites. IOV satellites show the largest STDs with 0.013 
cycles, while that of the two eccentric satellites is at the same level as other FOC satellites 
with 0.010 cycles. To further investigate the possible cause, more analyses including data 
usage rates, residual statistics are performed. 
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Figure 3 Time series of Galileo satellite WL FCBs for all the 17 satellites while E30 is selected as the refer- 
ence 
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Figure 4 STDs of daily satellite WL FCBs with 60 days. Averaged STDs are calculated for GPS and BDS, 
while an individual STD is calculated for each Galileo satellite. E30 is selected as the reference satellite for 


Galileo, while G01 and C09 are selected for GPS and BDS separately 
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Figure 5 HMW linear combinations of GPS, Galileo, and BDS dual frequency code measurements at sta- 


tion YEL2 on DoY 160, 2017 
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Figure 6 Usage rates of the float WL ambiguities. Averaged usage rates are calculated for GPS and BDS, 
while an individual usage rate is calculated for each Galileo satellite 


The float ambiguity observations, with residuals larger than 0.3 cycles, are rejected in 
FCB estimation. Figure 6 shows the usage rates of WL ambiguities for each satellite. We 
can see that almost all of the WL ambiguities are used for WL FCB estimation. The aver- 
aged usage rate is 96.2% for Galileo, while that for GPS and BDS were 95.3% and 90.1%, 
respectively. This indicates that the number of code multipath events for Galileo is 
smaller than for GPS, which may be ascribed to the advanced signal modulation of E5a. 
The usage rates of the two ECC satellites are comparable to other FOC satellites, while 
that for the IOV satellites are slightly lower. This indicates that the quality of FOC WL 
float ambiguities is slightly better than that of IOV satellites. The averaged data usage 
for BDS is about 90% which is much lower than that of Galileo and GPS. 
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Figure 7 Distribution of posterior residuals for Galileo (left), GPS (middle) and BDS (right) after WL FCB 


estimation 


The quality of FCB estimation can also be indicated by the posterior residuals of FCB 
estimation. Figure 7 shows the distribution of the residuals of WL FCB estimations. The 
mean values are -0.001, -0.003 and -0.003 cycles for Galileo, GPS, and BDS, respectively. 
The percentages of residuals smaller than 0.1 cycles are 86.5, 77.3 and 61.5% for Galileo, 
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GPS and BDS, respectively. The RMS of Galileo WL residuals is 0.071 cycles, while that 
for GPS and BDS are 0.089 and 0.117 cycles, respectively. The smaller and better distri- 
bution of Galileo residuals indicates a better quality of WL ambiguities. This further con- 
firms that the code multipath effects of Galileo are remarkably reduced. 


In order to further compare the WL FCB quality of IOV and FOC satellites, the RMS for 
each satellite in Galileo are presented in Figure 8. The figure also presents linear averages 
of RMS for each GPS and BDS satellite. It is found that the RMS of IOV satellites are 
significantly larger than other FOC satellites, while there is no significant difference be- 
tween FOC and ECC satellites. The average RMS of IOV satellites is 0.089 cycles, while 
that for FOC satellites is only 0.066 cycles. The percentage of residuals smaller than 0.1 
cycles is 76.4% and 88.9% for IOV and FOC satellites, respectively. These results indicate 
that the quality of FOC WL float ambiguities is better than that of IOV satellites. The 
worse performance of IOV satellites probably stems from the reduction of signal trans- 
mission power for IOV satellites. ESA imposed a reduction of 1.5 dB in the signal power 
of all four IOV satellites following a payload power problem of the fourth IOV (E20) in 
2014. The 1.5 dB power decrease results in an increase of approximately 15-20% of the 
thermal noise of the receiver, which roughly matches the observed increase in IOV WL 
FCB residuals. 
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Figure 8 RMS of WL residuals after WL FCB estimations. Averaged RMSs are calculated for GPS and BDS, 
while an individual RMS is calculated for each Galileo satellite 


From the above results, we can safely conclude that the overall performance of WL FCB 
can be expressed by the following inequality: 


FOC = ECC > IOV = GPS > BDS 


The WL FCB quality of the two ECC satellites is comparable to that of FOC satellites. 
Despite fewer observations acquired by receivers, the two satellites can be used for WL 
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ambiguity resolution. Concerning the lower data usage rates, larger residuals and larger 
STDs of WL FCB estimations for IOV satellites, the quality of IOV WL FCB is worse than 
that of FOC satellites. Despite the difference between the Galileo IOV and FOC satellites, 
the overall performance of Galileo WL FCB is better than that of GPS and BDS. 


4.2 Satellite NL FCB estimation 


Figure 9 presents the time series of Galileo/GPS/BDS satellite NL FCBs on DoY 160 in 
2017. For better visualization, the derived satellite FCBs are shifted with integer cycles. 
It can be seen that the derived satellite Galileo NL FCBs are only stable in tens of minutes. 
It is noted that the common mode shape of NL FCB variations is caused by the reference 
satellite as epoch-wise differencing is applied. 
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Figure 9 Time series of satellite NL FCBs for Galileo (upper), GPS (middle) and BDS (bottom), while satel- 
lite PRN E30, G10, C10 are selected as the references respectively. Each line represents one satellite, while 
the lines for reference satellites are zero and excluded. 


Figure 10 depicts the usage rate of float ambiguities for each Galileo satellite (average 
94.0%), as well as the average usage rates for GPS (97.6%) and BDS (80.4%). The usage 
rate of GPS NL float ambiguities is highest while that of BDS is worst. Since the weights 
of carrier phase measurements are 10000 times larger than that of code measurements, 
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the quality of the float ambiguities is dominated by the unmodeled errors in PPP. There- 
fore, the highest usage of GPS is reasonable since its precise product and models are 
currently state of the art. All Galileo satellites except E19 show similar usage rates rang- 
ing from 91.8 to 96.9%, indicating a similar quality level of NL ambiguities. The usage 
rate of E19 is only 83.8%, which is much smaller than that for other Galileo satellites. 
Inspired by Zaminpardaz and Teunissen (2017) who reported that for SEPTENTRIO re- 
ceivers the carrier-to-noise density ratio for E19 lies below the value of the other two IOV 
satellites for elevations higher than 60 degrees, we further calculated the carrier-to-noise 
density ratio, data usage rates for three receiver groups, namely LEICA (37 stations), 
SEPTENTRIO (41 stations) and TRIMBLE (83 stations). From the results in Figure 11, it 
can be seen that the carrier-to-noise density ratio of E19 lies below the other satellites for 
all the receiver groups. In addition, the data usage rates of E19 are lower than for the 
other satellites for all the three groups and no clear relationship with receiver type is 
found. Furthermore, we notice that the STD of satellite laser ranging (SLR) residuals for 
E19 (0.130 m) is larger than for E11 (0.092 m) and E12 (0.086 m) when performing SLR 
validation (Guo et al. 2017). This may indicate poor quality of precise satellite products 
for E19. Based on the above analysis, we suspect that the lower data usage rate of E19 is 
caused by the joint effect of poorer satellite products and larger carrier phase noise. 
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Figure 10 Usage rates of float NL ambiguities for satellite NL FCB estimation. Averaged usage rates are 


Percentage 


calculated for GPS and BDS, while an individual usage rate is calculated for each Galileo satellite 
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Figure 11 Carrier-to-noise density ratio C/NO of Galileo signals for three receiver groups, TRIMBLE (left), 
SEPTENTRIO (middle) and LEICA (right) 
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Figure 12 Distributions of posterior residuals for Galileo (left), GPS (middle) and BDS (right) after NL FCB 


estimation 


Figure 12 shows the distribution of the residuals of NL FCB estimations for the 60 days. 
The mean values are 0.000, 0.000 and -0.001 cycles for GPS, Galileo, and BDS, respec- 
tively. The percentages of residuals that are smaller than 0.1 cycles are 94.6, 91.4 and 
83.6% for GPS, Galileo, and BDS, respectively. The RMS of Galileo is 0.062, while that for 
GPS and BDS are 0.050 and 0.086 cycles, respectively. These results indicate that the 
quality of Galileo NL FCB estimations is slightly worse than that of GPS but better than 
that of BDS. As discussed above, the quality of NL FCB estimation could be affected by 
unmodeled errors. The lack of precise Galileo antenna PCO and PCV models both at 
satellite and user end certainly has influences on the PPP float ambiguities. The worst 
performance of BDS can be attributed to the poor quality of precise satellite products 
(Montenbruck et al. 2017), as well as the satellite-induced code multipath effects. 
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Figure 13 RMS of NL residuals after NL FCB estimations. Averaged RMSs are calculated for GPS and BDS, 
while an individual RMS is calculated for each Galileo satellite 


In order to compare the NL FCB quality of IOV and FOC satellites, the RMSs for each 
Galileo satellite are presented in Figure 13. It can be found that there is no significant 
difference among FOC, IOV and ECC satellites. The average RMSs are 0.061, 0.065 and 
0.058 cycles for FOC, IOV, and ECC, respectively. The percentages of residuals that are 
smaller than 0.1 cycles are 91.9, 90.3 and 92.7% for FOC, IOV, and ECC, respectively. 
These results indicate that the satellites in the incorrect orbital planes behave not worse 
than those in nominal FOC and IOV orbits in terms of PPP. We can further deduce that 
the orbit and clock quality of ECC satellites from MGEX precise product are at a similar 
level as that of FOC and IOV. The largest RMS is found for satellite E19 with values at 
about 0.069 cycles. Concerning the lower data usage of this satellite, the quality of its NL 
FCB is the worst. From the results of the usage rates of NL float ambiguities and distri- 
butions of the residuals, we conclude that the overall performance of Galileo NL FCB is 
worse than that of GPS but better than BDS. It is noted that the ranking is temporary, 
and promising improvements can be expected for Galileo and BDS considering their 
rapid developments. 


4.3 PPP AR solution 


In order to validate our FCB estimates, as well as to assess the performance of PPP AR, 
a subset of MGEX network stations is processed in PPP AR mode. These test stations 
have been excluded from the FCB estimation. Due to the incomplete constellation of Gal- 
ileo, the number of observable satellites differs with respect to geographical locations. 20 
MGEX stations, as shown in Figure 14, with relatively higher percentage of Galileo ob- 
servations are selected. GPS and Galileo observations from DoY 155 to 164 are processed. 
For static solutions, the 24 hours observations are divided into eight three-hour sessions. 
For each session, the epochs with less than 4 visible Galileo satellites are deleted. Then 
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the sessions with data integrity less than 80 percent are removed. After the pre-pro- 
cessing, there are 1484 sessions. For pseudo-kinematic solutions, the 24 hours observa- 
tions are regarded as one session. Galileo-only solutions, GPS-only solutions, as well as 
GPS/Galileo-combined solutions, are calculated for each session. The positional biases 
of static Galileo PPP float solutions and AR solutions on DoY 160, 2017 are presented in 
Figure 15. The statistics of all sessions in the ten days are provided in Table 2. 


Figure 14 Geographical distribution of the selected 20 MGEX stations 
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Figure 15 Absolute NEU biases of static Galileo PPP float and AR solutions with three-hour sessions on 
DoY 160, 2017. In total, there are 98 sessions for the 20 MGEX stations 
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Table 2 Statistics of static PPP float and AR solutions with three-hour sessions for 20 stations in ten days 


[Unit: m] 

Solut. No. of successful sessions East North Up 
GAL-Float 994/1484 0.021 0.011 0.027 
GAL-AR 994/1484 0.007 0.006 0.021 
GPS-Float 1422/1484 0.015 0.007 0.020 
GPS-AR 1422/1484 0.006 0.005 0.018 
GAL/GPS-Float 1442/1484 0.012 0.007 0.018 
GAL/GPS-AR 1442/1484 0.006 0.005 0.017 


It can be seen that the accuracy of the Galileo-only PPP float solution is several centime- 
ters with the current constellation. After the ambiguity resolutions, the accuracy is sig- 
nificantly improved. The most significant improvement is found for the east component, 
from 2.1 to 0.7 cm, corresponding to an enhancement of 67%. For GPS, the improvement 
of the north component by PPP AR is insignificant as the accuracy of the north compo- 
nent in the float solution is relatively high. However, in Galileo, we also observe an ob- 
vious improvement for the north component, from 1.1 to 0.6 cm with an enhancement of 
45%. This is because the accuracy of the north component in Galileo float solution is low 
due to bad geometry. The improvement for the vertical component is about 22%, which 
is less significant than for the horizontal components. The possible reason is that the 
vertical component is tightly coupled with receiver clock and zenith wet delay parame- 
ters, which are simultaneously estimated. 


Compared with the GPS-only solutions, both float and AR solutions of Galileo are 
slightly worse, although it has been proven that the quality of Galileo WL FCB is better 
than GPS. This can partly be ascribed to the smaller number of available satellites. It can 
also be ascribed to the worse quality of Galileo NL FCB as the accuracy of NL FCB di- 
rectly affects the accuracy of PPP AR. 
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Figure 16 NEU biases of Galileo, GPS, GPS+Galileo PPP AR in static mode, at station YEL2 on DoY 160, 
2017 


Figure 16 shows a typical example for standalone single- and combined PPP AR solu- 
tions in static mode. One can see that the Galileo-only solution requires more time to 
achieve a reliable accuracy, while the convergence time of GPS-only solutions is much 
shorter. The fastest convergence can be achieved in GPS/Galileo PPP AR. Although the 
final accuracy of the GPS/Galileo combined solution is almost the same as that of the 
GPS-only solution, integrating Galileo with GPS increases the fixing rate and results in 
a more stable solution. This is further confirmed by the reduction of incorrectly fixed 
sessions (see Table 2), where the success rate increases from 95.8% to 97.2% when adding 
Galileo to GPS. 
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Figure 17 NEU biases of kinematic Galileo PPP float and AR solutions at station YEL2 on DoY 160, 2017 


Figure 17 presents the kinematic positional biases for station YEL2 on DoY 160, 2017. 
Epochs with less than five valid satellites have been deleted. It can be seen that the accu- 
racy of all three components is significantly improved by ambiguity resolution. In addi- 
tion, it is observed that the convergence time is also shortened. And yet the time to first 
fix (TTFF) is still larger than two hours as a result of the limited number of observable 
satellites. 
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Figure 18 NEU biases of kinematic GPS, and GPS/Galileo PPP AR at station YEL2 on DoY 160, 2017 


Although Galileo-only solutions are limited by the number of observable satellites, add- 
ing Galileo to GPS can enhance the performance of kinematic PPP solutions. Figure 18 
shows a typical example for GPS-only and GPS/Galileo combined PPP AR solutions. As 
we can see, integrating Galileo with GPS increases the accuracy of kinematic AR solu- 
tions, especially when the accuracy of the GPS-only solution is bad. The combined solu- 
tion can achieve an accuracy of 1-2 cm in horizontal and 4-6 cm in vertical components. 


Table 3 Statistics of kinematic PPP float and AR solutions for 20 stations in ten days [Unit: m] 


16 


Solut. East North Up 
GAL-Float 0.027 0.022 0.050 
GAL-AR 0.011 0.012 0.036 
GPS-Float 0.018 0.013 0.035 
GPS-AR 0.009 0.009 0.030 
GAL/GPS-Float 0.015 0.011 0.030 
GAL/GPS-AR 0.008 0.008 0.028 


Pseudo-kinematic PPP AR solutions are conducted for each station in the ten days. The 
statistics of positional biases are presented in Table 3. For Galileo-only solutions, the 
RMSs of the east, north and up components are improved by 59%, 45%, and 28% respec- 
tively, while for GPS-only solutions the enhancements are 50%, 31% and 14%. Similar to 
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the static results, the improvements for Galileo are more significant than for GPS. Dual- 
system PPP AR can further improve the kinematic positional accuracy. 


5 Conclusions 


The feasibility of Galileo PPP ambiguity resolution with the current constellation has 
been demonstrated. The satellite WL/NL FCBs are estimated from globally distributed 
MGExX stations. Results indicate that the quality of Galileo WL FCB is better than for GPS 
and BDS in terms of data usage rate, residual distribution, as well as STD of daily esti- 
mates. We attribute the good quality of Galileo WL FCB to its advanced signal modula- 
tion, AltBOC, which significantly compresses the multipath effect for code measure- 
ments. Within the Galileo constellation, the quality of FOC WL FCB is much better than 
for IOV satellites. The poorer performance of IOV satellites WL FCB is a result of a re- 
duction in the satellite transmit signal power. The performance of the two satellites with 
highly eccentric orbits is comparable to other FOC satellites but having a smaller number 
of observations. As for NL FCB, the quality of Galileo NL FCB is slightly worse than that 
of GPS but better than that of BDS. Since the accuracy of NL FCB estimation is dominated 
by unmodeled errors in float PPP, the worse quality of Galileo NL FCB is likely caused 
by the imperfect PCO and PCV models used in this study. Within Galileo, the NL FCB 
quality of FOC and IOV (except E19), as well as the two eccentric satellites, shows no 
significant difference in terms of data usage rates and residuals. The reason for the worse 
performance of E19 is still not clear. On the one hand, it cannot, or at least not fully, be 
ascribed to the signal transmission power as the power difference between FOC and IOV 
E11/E12 is larger than the difference between E11/E12 and E19. On the other hand, the 
worse quality of E19 satellite orbit, indicated by SLR residuals, could also degrade NL 
FCB estimations. This issue remains an open question and deserves further investiga- 
tion. For static Galileo PPP AR with three-hour sessions, the positional biases can be re- 
duced by 67, 45 and 22% for east, north and up components, respectively. In addition, 
PPP AR also improves the kinematic solutions by 59, 45 and 28% for east, north and 
vertical components, respectively. Although the performance of Galileo PPP AR is still 
worse than that of GPS, a promising improvement can be expected in the near future as 
the Galileo FOC satellite metadata, including satellite mass, attitude law, PCO and PCV, 
are published on 2017-10-06 (https://www.gsc-europa.eu/support-to-developers/galileo- 
satellite-metadata). This informa-tion is expected to improve the precise orbit determi- 
nation and PPP solution. 


Acknowledgments Thanks go to IGS-MGEX and GFZ for providing GNSS data and products. Thanks also 
go to Maorong GE from GFZ for valuable discussions. Guorui Xiao is supported by the China Scholarship 
Council. This work was supported in part by the National Natural Science Foundation of China (grants nos. 
41674016, 41274016 and 41604024). Finally, we highly acknowledge the comments by two anonymous re- 
viewers which helped to improve the submitted version of the paper. 


81 


Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution 


References 


21. 


82 


Collins P, Lahaye F, Héroux P, Bisnath S (2008) Precise point positioning with ambiguity resolution 
using the decoupled clock model. In: Proc. ION GNSS 2008, Institute of Navigation, Savannah, Georgia, 
September 16-19 2008. pp 1315-1322 

Ge M, Gendt G, Rothacher Ma, Shi C, Liu J (2008) Resolution of GPS carrier-phase ambiguities in precise 
point positioning (PPP) with daily observations. Journal of Geodesy 82(7):389-399 

Geng J, Meng X, Dodson AH, Teferle FN (2010) Integer ambiguity resolution in precise point position- 
ing: method comparison. Journal of Geodesy 84(9):569-581 

Geng J, Shi C (2016) Rapid initialization of real-time PPP by resolving undifferenced GPS and 
GLONASS ambiguities simultaneously. Journal of Geodesy 91(4):361-374 

Gu S, Lou Y, Shi C, Liu J (2015) BeiDou phase bias estimation and its application in precise point posi- 
tioning with triple-frequency observable. Journal of Geodesy 89(10):979-979 

Guo F, Li X, Zhang X, Wang J (2017) Assessment of precise orbit and clock products for Galileo, BeiDou, 
and QZSS from IGS Multi-GNSS Experiment (MGEX). GPS Solutions 21(1):279-290 

Hatch R (1982) The synergism of GPS code and carrier measurements. In: Proceedings of the third in- 
ternational symposium on satellite doppler positioning at Physical Sciences Laboratory of New Mexico 
State University, February 8-12, pp 1213-1231 

Laurichesse D, Mercier F, Berthias J-P, Broca P, Cerri L (2009) Integer ambiguity resolution on undiffer- 
enced GPS phase measurements and its application to PPP and satellite precise orbit determination. 
Navigation 56(2):135-149 

Li P, Zhang X (2015) Precise point positioning with partial ambiguity fixing. Sensors 15(6):13627-13643 


. LiP, Zhang X, Ren X, Zuo X, Pan Y (2015) Generating GPS satellite fractional cycle bias for ambiguity- 


fixed precise point positioning. GPS Solutions 20(4):771-782 


. Li P, Zhang X, Guo F (2017a) Ambiguity resolved precise point positioning with GPS and BeiDou. Jour- 


nal of Geodesy 91(1):25-40 


. Li X, Li X, Yuan Y, Zhang K, Zhang X, Wickert J (2017b) Multi-GNSS phase delay estimation and PPP 


ambiguity resolution: GPS, BDS, GLONASS, Galileo. Journal of Geodesy 92(6):579-608 


. Li X, Zhang X (2012) Improving the estimation of uncalibrated fractional phase offsets for PPP ambigu- 


ity resolution. Journal of Navigation 65(3):513-529 


. Liu Y, YeS, Song W, Lou Y, Gu S (2017) Rapid PPP ambiguity resolution using GPS+ GLONASS obser- 


vations. Journal of Geodesy 91(4):441-455 


. Loyer S, Perosanz F, Mercier F, Capdeville H, Marty J-C (2012) Zero-difference GPS ambiguity resolu- 


tion at CNES-CLS IGS Analysis Center. Journal of Geodesy 86(11):991-1003 

Melbourne WG (1985) The case for ranging in GPS-based geodetic systems. In: Proceedings of the first 
international symposium on precise positioning with the Global Positioning System, Rockville, 1985. 
pp 373-386 

Montenbruck O et al. (2017) The multi-GNSS experiment (MGEX) of the international GNSS service 
(IGS) — Achievements, prospects and challenges. Advances in Space Research 59(7):1671-1697 


. Prange L, Orliac E, Dach R, Arnold D, Beutler G, Schaer S, Jaggi A (2017) CODE 's five-system orbit and 


clock solution —the challenges of multi-GNSS data analysis. Journal of Geodesy 91(4):345-360 


. Shi J, Gao Y (2014) A comparison of three PPP integer ambiguity resolution methods. GPS Solutions 


18(4):519-528 


. Sonica K, Prange L, Kamierski K, Bury G, Drodewski M, Zajdel R, Hadas T (2017) Validation of Galileo 


orbits using SLR with a focus on satellites launched into incorrect orbital planes. Journal of Geodesy 
92(2):131-148 

Takasu T, Yasuda (2010) A Kalman-filter-based integer ambiguity resolution strategy for long-baseline 
RTK with ionosphere and troposphere estimation. In: ION NTM 2010. pp 161-171 


Chapter 3 Estimating and assessing Galileo satellite fractional cycle bias for PPP ambiguity resolution 


22. Tegedor J, Liu X, Jong K, Goode M, Ovstedal O, Vigen E (2016) Estimation of Galileo Uncalibrated 
Hardware Delays for Ambiguity-Fixed Precise Point Positioning. Navigation 63(2):173-179 

23. Teunissen PJG, Jonge PJ, Tiberius CCJM (1997) Performance of the LAMBDA method for fast GPS am- 
biguity resolution. Navigation 44(3):373-383 

24. Teunissen PJG, Joosten P, Tiberius CCJM (1997) Geometry-free ambiguity success rates in case of partial 
fixing. In: ION NTM 1999. pp 25-27 

25. Teunissen PJG, Khodabandeh A (2015) Review and principles of PPP-RTK methods. Journal of Geodesy 
89(3):217-240 

26. Uhlemann M, Gendt G, Ramatschi M, Deng Z (2016) GFZ Global Multi-GNSS Network and Data Pro- 
cessing Results. Springer International Publishing 

27. Wang M, Chai H, Li Y (2017) Performance analysis of BDS/GPS precise point positioning with undiffer- 
enced ambiguity resolution. Advances in Space Research 60(12):2581-2595 

28. Wanninger L, Beer S (2015) BeiDou satellite-induced code pseudorange variations: diagnosis and ther- 
apy. GPS Solutions 19(4):639-648 

29. Wiibbena G (1985) Software developments for geodetic positioning with GPS using TI-4100 code and 
carrier measurements. In: Proceedings of the first international symposium on precise positioning with 
the global positioning system, Rockville, 1985. pp 403-412 

30. Xiao G, Sui L, Heck B, Zeng T, Tian Y (2018) Estimating satellite phase fractional cycle biases based on 
Kalman filter. GPS Solutions 22(3):82. https://doi.org/10.1007/s10291-018-0749-3 

31. Yi Wet al. (2017) Improved method to estimate undifferenced satellite fractional cycle biases using net- 
work observations to support PPP ambiguity resolution. GPS Solutions 21(3):1369-1378 

32. Zaminpardaz S, Teunissen PJG (2017) Analysis of Galileo IOV+ FOC signals and E5 RTK performance. 
GPS Solutions 21(4):1855-1870 


Author Biographies 


Guorui Xiao is currently a Ph.D. candidate at Geodetic Institute, Karlsruhe Institute of Technol- 
ogy (KIT). He received his B.Sc. degree from the School of Geodesy and Geomatics in Wuhan 
University in 2011. His current research focuses on multi-frequency and multi-constellation 
GNSS precise positioning and applications. 


Pan Li is a postdoctoral researcher at GFZ, Germany. He obtained his Ph.D. degree in 2016 from 
School of Geodesy and Geomatics, Wuhan University, China. His current research focuses mainly 
on GNSS precise point positioning and ambiguity resolution. 


Lifen Sui is a professor at Zhengzhou Institute of Surveying and Mapping. She obtained her 
Ph.D. degree from Zhengzhou Institute of Surveying and Mapping in 2001. Her main research 
interests include surveying data processing and GNSS positioning. 


Bernhard Heck received the Ph.D. degree and the Habilitation in Geodesy from the University 
of Karlsruhe (now: KIT), Germany, in 1979 and 1984, respectively. Since 1991 he fills the position 
of Professor and Chair of Physical and Satellite Geodesy at KIT. His research interests include 
geometrical geodesy, deformation analysis, gravity field modelling, and GNSS positioning. 


Harald Schuh obtained his Ph.D. from the University of Bonn, Germany and is currently the 
Director of the Department 1 “Geodesy”, German Research Centre for Geosciences (GFZ). He is 


83 


Multi-frequency and multi-GNSS PPP phase bias estimation and ambiguity resolution 


also the President of the International Association of Geodesy (LAG) since 2015. His research in- 
terests are the space geodetic techniques in particular Very Long Baseline Interferometry and 
GNSS and their applications. 


84 


Chapter 4 

A Unified Model for Multi- 
Frequency PPP Ambiguity 
Resolution 


Guorui Xiao’* - Pan Li**- Yang Gaot- Bernhard Heck! 


! Geodetic Institute, Karlsruhe Institute of Technology, D-76131 Karlsruhe, Germany; bern- 
hard.heck@kit.edu 

? Zhengzhou Institute of Surveying and Mapping, Zhengzhou 450052, China 

3 The German Research Centre for Geosciences (GFZ), D-14473 Potsdam, Germany 

4 Department of Geomatics Engineering, University of Calgary, Calgary, AB T2N 1N4, Canada; ygao@ucal- 
gary.ca 

* Correspondence: xgr@whu.edu.cn (G.X.); panli@gfz-potsdam.de (P.L.) 


Submitted: 11 December 2018; Accepted: 7 January 2019; Published: 10 January 2019 


Remote Sensing 11(2):116, 2019. DOI:10.3390/rs11020116 


This is an author-created version of the article with permission from MDPI. The final publication is availa- 
ble at www.mdpi.com 


Abstract: With the modernization of Global Navigation Satellite System (GNSS), triple- or multi- 
frequency signals have become available from more and more GNSS satellites. The additional 
signals are expected to enhance the performance of precise point positioning (PPP) with ambigu- 
ity resolution (AR). To deal with the additional signals, we propose a unified modeling strategy 
for multi-frequency PPP AR based on raw uncombined observations. Based on the unified model, 
the fractional cycle biases (FCBs) generated from multi-frequency observations can be flexibly 
used, such as for dual- or triple- frequency PPP AR. Its efficiency is verified with Galileo and 
BeiDou triple-frequency observations collected from globally distributed MGEX stations. The es- 
timated FCB are assessed with respect to residual distributions and standard deviations. The ob- 
tained results indicate good consistency between the input float ambiguities and the generated 
FCBs. To assess the performance of the triple-frequency PPP AR, 11 days of MGEX data are pro- 
cessed in three-hour sessions. The positional biases in the ambiguity-fixed solutions are signifi- 
cantly reduced compared with the float solutions. The improvements are 49.2%, 38.3%, and 
29.6%, respectively, in east/north/up components for positioning with BDS, while the correspond- 
ing improvements are 60.0%, 29.0%, and 21.1% for positioning with Galileo. These results confirm 
the efficiency of the proposed approach, and that the triple-frequency PPP AR can bring an obvi- 
ous benefit to the ambiguity-float PPP solution. 


Keywords: Galileo - BeiDou - Precise point positioning - Integer ambiguity resolution - Triple- 
frequency - Fractional cycle bias 
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1 Introduction 


Precise point positioning (PPP) has found increased applications due to its cost-effec- 
tiveness, global coverage, and high accuracy (Kouba and Héroux 2001; Zumberge et al. 
1997). Usually PPP is able to achieve a positional accuracy of 10 cm after a convergence 
time of 30 min (Cai et al. 2015). The integer ambiguity resolution (AR) technique is ex- 
pected to further enhance the accuracy and shorten the convergence time (Gabor and 
Nerem 1999; Gao and Shen 2002). In addition, (Katsigianni et al. 2018) shows that the 
Galileo orbit determination could be improved when employing AR in multiple Global 
Navigation Satellite System (GNSS) data processing. However, the uncalibrated phase 
delays (UPDs) originating from satellites and receivers destroy the integer nature of PPP 
ambiguities. By determining the UPDs, to be estimated as fractional cycle biases (FCBs) 
at the server end and applying them at the user end, the PPP integer ambiguity resolu- 
tion could become feasible (Ge et al. 2008; Li et al. 2015a; Xiao et al. 2018b). Similarly, the 
decoupled clock model (Collins et al. 2008) and the integer phase clock model 
(Laurichesse et al. 2009) were developed. These PPP AR techniques have been proven 
equivalent in theory (Shi and Gao 2014; Teunissen and Khodabandeh 2015), and the po- 
sitional biases have been demonstrated to be minimal (Geng et al. 2010). Beside GPS, 
PPP AR has been extended to GLONASS (Geng and Shi 2016; Liu et al. 2017b; Yi et al. 
2017), BeiDou Navigation Satellite System (BDS) (Li et al. 2017a; Liu et al. 2017a), Galileo 
(Tegedor et al. 2016; Xiao et al. 2019b), and multi-GNSS (Li et al. 2017b). 


As to the functional models used for PPP, the dual-frequency ionospheric-free (IF) com- 
bination is routinely employed (e.g., in the above-mentioned research). However, with 
emerging BDS and Galileo, as well as the modernization of GPS and GLONASS, various 
types of multi-frequency observables become available (Montenbruck et al. 2017). The 
choice of optimum combinations then becomes practically difficult given the diversity 
of equipment (Schönemann et al. 2011). In addition, the IF combination will amplify the 
measurement noise level by a factor of about 3, which will degrade the performance of 
the position solution. As a result, the PPP model based on uncombined measurements, 
in which the individual signal of each frequency is treated as an independent observable, 
has drawn increasing interest in the GNSS community (Gu et al. 2015b). Its efficiency 
has already been confirmed in terms of convergence time and precision for single-fre- 
quency PPP (Lou et al. 2016), multi-GNSS PPP (Chen et al. 2015), as well as PPP-RTK 
(Feng et al. 2013b; Odijk et al. 2016b). Moreover, this approach has been tested effectively 
for ionospheric modelling (Tu et al. 2013; Zhang et al. 2012), differential code bias (DCB) 
estimation (Liu et al. 2018; Shi et al. 2015), and low earth satellite orbit determination 
(Zehentner and Mayer-Gürr 2016). 


Compared to the well-developed IF PPP model, the uncombined PPP model, called un- 
combined PPP in the sequel, requires more investigation, especially in the case of multi- 
frequency processing and ambiguity resolution. First, how to model and constrain the 
ionospheric delay has a crucial impact on performance using the uncombined PPP. For 
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example, it has been demonstrated that a white noise model is not adequate to capture 
the characteristics of the ionospheric delay. The external constraints developed from the 
ionospheric products, such as the IGS global ionosphere maps, are also not accurate 
enough to completely separate the ionospheric effects from the ambiguity parameters 
(Gu et al. 2015b). The influence of the ionospheric effects on the ambiguity fixing there- 
fore must be reduced. Second, the method to deal with the DCB errors is more problem- 
atic with the uncombined PPP (Guo et al. 2015) than the dual-frequency IF PPP, since 
the latter can cancel out the DCB biases. The problem of partial assimilation of the code 
bias (DCB) into phase bias (FCB) should also be carefully considered. Third, the uncombined 
PPP approach was proposed to deal with multi-GNSS and multi-frequency signals, so a 
generalized FCB estimation and AR method (Li et al. 2018), which is extendable to dual-, 
triple-, and multi- frequency, should be proposed. 


Li et al. (2013a) verified the feasibility of the uncombined PPP AR with refined iono- 
spheric models. The ionospheric delay was constrained from a priori spatial-temporal 
information and ionospheric products. The GPS dual-frequency ambiguities were fixed 
sequentially in the forms of wide-lane (WL)/narrow-lane (NL), which followed the con- 
vention of IF PPP AR. Gu et al. (2015a) further testified the uncombined PPP AR with 
BDS triple-frequency observations. The extra-wide-lane (EWL) and WL ambiguities 
were successfully fixed, whereas the B1 ambiguities were kept as float values. In addi- 
tion, the performance was further limited by the satellite-induced multipath effects 
(Wanninger and Beer 2015). Li et al. (2018) proposed a unified FCB estimation and PPP 
AR method, which is extendable to multi-frequency uncombined PPP. The FCBs on each 
frequency were directly estimated from the raw float ambiguities derived from triple 
frequency observables. The model showed a great potential for multi-frequency uncom- 
bined PPP AR, although its DCB strategy may not be optimal. The satellite DCBs, to- 
gether with the receiver DCB, were estimated as unknowns, and as a result the number 
of unknown parameters was increased. Given that the satellite GNSS DCB product is 
currently available on a routine basis (Wang et al. 2016), it would be beneficial to make 
use of these products. In addition, validating the method with Galileo observations is of 
interest considering the recent and rapid development of Galileo. 


The aim of this study is to develop a unified modeling strategy for multi-frequency and 
multi-GNSS uncombined PPP AR. The unified model is able to generate consistent FCB 
products and perform PPP AR for multi-frequency PPP. The proposed approach will be 
first described, and its effectiveness will be verified with Galileo and BDS triple-fre- 
quency observations collected from the globally distributed Multi-GNSS EXperiment 
(MGEX) stations. The estimated FCB are assessed with respect to residual distributions 
and standard deviations, followed by an evaluation of the performance improvements 
in Galileo and BDS triple-frequency PPP AR. Finally, the results are summarized, and 
an outlook for future research is presented. 
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2 Methodology 


The proposed uncombined PPP mode will be first described in this section, followed by 
a description of our FCB estimation strategy. With the obtained FCB products, the un- 
combined PPP AR algorithm at the user end is then elaborated. 


2.1 Uncombined PPP Float Ambiguity Model 


In the classic GNSS dual-frequency PPP, the first-order ionospheric delay is eliminated 
by the formation of the IF combination (Zumberge et al. 1997). In the uncombined PPP 
model, the ionospheric delay is directly estimated. For a satellite s observed by receiver 
r, the corresponding raw pseudo-range and carrier phase observation equations can be 
expressed as (Leick et al. 2015) 


ee e 
(1) 


nf = pe + dt, — dt’ + dT — ap: dlp, Aj(NZ + By — Bj) + £o, 


where the subscript f = (1,2,3,---) refers to a specific carrier frequency, superscript s re- 
fers to a specific satellite; p? indicates the geometric distance between the satellite and 
receiver; dt, and dt? are the clock errors of receiver and satellite; dT is the slant tropo- 
spheric delay; dl; is the slant ionospheric delay on the first carrier frequency and ay = 
Àj / À is the carrier frequency-dependent factor; D, p and Dj are the receiver and satel- 
lite specific code hardware delays; Ay and N; , are the wavelength in meter and integer 
ambiguity in cycle; B, and Bj are the receiver-dependent and satellite-dependent un- 
calibrated phase delays; £p, and £5, are the pseudo-range and carrier phase measure- 
ment noise, respectively. Note that the higher-order ionospheric effects are neglected, as 
they have limited influence on the performance of ambiguity resolution (Hadas et al. 
2017). 

Another important difference between the IF PPP and the uncombined PPP model is the 
strategy to deal with the DCB. The DCB is not of concern in IF PPP as the IF combination 
is also used for precise clock generation, which implies that the DCB could be fully ab- 
sorbed by other parameters or simply cancelled out in the IF PPP (Dach et al. 2015). But 
this is not the case in the uncombined PPP, especially with multi-frequency observations. 
Conventionally, precise orbit and clock products from the IGS analysis center are used 


to remove satellite orbit and clock errors. During the generation of precise clock prod- 

a2 E 
a,-1 4 
the clock offset dt? in accordance with the IGS analysis convention. After applying the 
GNSS precise satellite clock products, Equation (1) can be rewritten as 


ucts, the pseudo-range IF hardware delay bias Dj; = u a D3 is assimilated into 
T 


5, = pj dt, dtre dT + af: dls, + Doy — Dj + Di £p, 
(2) 


np = Pp + dt, — dtre + dT — ap: dl? + Ap (Nop + By — Bj) + Die £o, 
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This linear system is rank-deficient due to the DCB parameters. For dual-frequency un- 
combined PPP processing, the singularities can be eliminated by a re-parameterization 
process. This is accomplished based on the fact that the DCB parameters can be sepa- 
rated into satellite-related, receiver-related, and frequency-related parts, and therefore 
can be fully absorbed, respectively, by satellite clock, receiver clock, and ionospheric pa- 
rameters (Zhang et al. 2012). This method is efficient for dual-frequency processing but 
becomes complicated when facing multi-frequency observations. It is also possible to 
estimate these DCB parameters in advance. This is typically done by employing a net- 
work of receivers and imposing a zero-mean constraint. This option is complicated and 
not suitable for single receivers. In our study, we propose estimation of the receiver DCB 
and correction of the satellite DCB with existing multi-GNSS DCB products (Wang et al. 
2016). Taking triple-frequency observations as an example, the correction equation can 
be deduced as (Guo et al. 2015) 


J 1 
Pra = pr + dt, — dtre +dT + a: dl, umi DCB, + Ep, 
2— 


_ a 
PS, = pj + dt, — dt$,, + dT + a; di}, — DCBi + DCB,.12 + £p, (3) 
= 
u 1 
Fa = p? + dt, — dtöre + dT + az dl}, — DCB}; — —— DCB}, + DCByas + Ep, 
"T 
where DCB$, = D$ — Dj, DCB,12 = D; — D,4, and dt, = dt, + D,. The DCBS, and 
DCB}; can be obtained from multi-GNSS DCB products, while DCB, 4; and DCB,.45 are 


estimated as daily constant parameters. Similarly, the phase equations can be rewritten 
as 


rf = ps + dt, — thre + dT — ag: dlz, + Ay Nr p E Ty 


where the ambiguity can be re-parameterized as 


Dra 
b, ¢ = Beg —— 
er (5) 
S S iF 
om 
and the estimable parameters are 
= S —S —S —s 
x= |x y z dt. dT I Dri2 Dras Nra Ny2 LAM (6) 


Compared with the model in Li et al. (2018), the structure of the unknown parameters, 
except (x, y, Z), is different, due to the different strategies of DCB correction. In Li et al. 
(2018), the dual-frequency DCBs are absorbed by other parameters, whereas the third 
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frequency DCBs are estimated. In this study, the satellite DCBs for all the three frequen- 
cies are corrected with existing DCB products (Wang et al. 2016) and the receiver DCBs 
are estimated. Consequently, our ionospheric parameters will not be biased by DCBs, 
which is beneficial for ionospheric modelling. In addition, for single stations with n ob- 
servable satellites, the number of DCB parameters to be estimated in their model is n, 
while it is 2 in our model. The degree of freedom of our model is larger, which could 
increase the redundancy and robustness of the positioning solutions. The estimated am- 
biguity parameter is a combination of the integer ambiguity, the corresponding code 
hardware delays, and the uncalibrated carrier phase delays at both receiver and satellite 
ends. In order to recover its integer property, these biases, i.e., satellite FCB b? and re- 
ceiver FCB b, f, must be accounted for. Normally, the receiver FCB is not of concern as 
it can be eliminated when performing single differences of observations between satel- 
lites. The satellite FCB, however, must be estimated at the server end and broadcasted to 
the users. 


2.2 FCB Estimation Strategy 


In dual-frequency IF PPP, the float ambiguity is usually decomposed into WL/NL forms 
in order to recover the integer property (Ge et al. 2008). This is partly because the IF 
combination of L1/L2 ambiguities is, in essence, not an integer. Another reason is that 
the WL ambiguities possess a relatively longer wavelength and are less correlated, there- 
fore can be easily fixed. For uncombined PPP AR, it is also important to form combina- 
tions of raw ambiguities. On the one hand, the estimated raw float ambiguities are 
strongly correlated. On the other hand, the raw float ambiguities are quite sensitive to 
unmodeled ionospheric errors (Gu et al. 2015b). Therefore, the combinations with longer 
wavelengths and lower ionospheric delays are preferred. The coefficients must be inte- 
gers in order to preserve the integer nature of ambiguities. In addition, these combina- 
tions should be independent to avoid rank-deficiency. While the Least-squares AMBi- 
guity Decorrelation Adjustment (LAMBDA) method can be used to automatically search 
for the optimal linear combinations of ambiguities (Li et al. 2018; Teunissen et al. 1997), 
the classic extra-/wide-lane ambiguities (EWL/WL) were found to perform equally well 
and were used in our work to simplify the algorithm. For triple-frequency observations, 
it is easy to find two optimal combinations, e.g., one EWL and one WL combination or 
two WL combinations. The searching of the third combination, however, is much more 
difficult. From the systematic investigation of triple-frequency combinations, it is found 
that (4, -3, 0) is a good compromise between ionospheric reduction and noise amplifica- 
tion (Li et al. 2018). Concerning the properties of the above combinations shown in Table 
1, they are denoted as NL/WL/EWL without specific explanation herein. In addition, 
they are used for both BDS and Galileo triple-frequency observations for simplicity 
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S —S 


A yo 1 4 —3 0 " n 

Nyic2|= |1 -1 0||N;.2 (7) 
—s 1 0 —1l|—s 

N r,LC3 N ra 


Substituting (5) into the above system produces the basic model for estimating FCBs. 
Since they have the same structure, a general expression can be formulated as 


—s oe 
S <2 S S 
rLC = Nyc D Nr ic = drc Ze dic (8) 


for all linear combinations, N, tc denotes the float combined ambiguities; Ñ; denotes the 
integer part of Ny ic; dr ıc and dj. denotes the receiver and satellite FCBs; Rz; repre- 
sents the FCB measurements. For each linear combination, a set of equations in the form 
of (8) can be generated, based on a network of reference stations. Suppose that there are 
n satellites tracked by m reference stations, the system of equations can be expressed as 


Rlıc Nic — Nine 15 0 -1 OJrd 

: : AE 1,LC 

y —n = : : : . : : 

Ric Nic — Nic 1-00 - -1 4 

: k : : . . . i m,LC 

"wm ie ° ve Xe t epo € Xo (9) 
N X d =f iw LC 

Rac] |Wmic-Whec] | voi Cb ov ON 

i : . . . . . = d?. 

Rm,ıc AD o Ow 1 QO = -1 


n 
Nanc m Nm 


where LC stand for the linear combinations (LC1, LC2, LC3, ---). The obtained system is 
singular on both sides of the equations. For the left side, N,. can be determined by 
rounding N;,c, assuming that the float ambiguities are precisely estimated. For the right 
side, one arbitrarily combined FCB should be set to zero. For all the linear combinations, 
we always set the combined FCB of the last satellite to zero, i.e., G32/E30/C14. Note that 
the FCB measurements R} ¿cç from different stations may differ with +1 cycle. This is due 
to the rounding process and can be adjusted with the strategy described in Xiao et al. 
(2018b). In this way, the system of equations can be solved. With the obtained combined 
FCB, we are able to calculate the FCB of the raw L1/L2/L3 carrier frequency 


d| [4-s3 oT!ldaı 
dij|2|1 -1 0| Id (10) 
dj ı 0 -1l [des 


The transformation from combined FCBs to raw FCBs is important, as it provides more 
flexibility to the users. With the raw FCBs, users are able to choose their own linear com- 
binations of observations, formulate the corresponding combined FCB, and conduct PPP 
AR. This representation allows interoperability if the server and user sides implement 
different AR methods. In addition, the raw FCB is suitable for the State Space Represen- 
tation (SSR) of Radio Technical Commission for Maritime services (RTCM) (Weber et al. 
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2005), where one phase bias per phase observable is broadcasted instead of making spe- 
cific combinations. 


Table 1 Properties of GPS, Galileo, and BDS triple-frequency linear combinations. 


GNSS Coefficients Wavelength [meter] Ionospheric Delay [cycle] Noise [cycle] 
(4, -3, 0) 0.114 0.150 5.0 
GPS (1, -1, 0) 0.862 -0.283 1.414 
(1, 0, -1) 0.751 -0.339 1414 
(4, -3, 0) 0.108 -0.017 5.0 
Galileo (1,-1,0) 0.751 -0.339 1.414 
(1,0,-1) 0.814 -0.305 1.414 
(4, -3, 0) 0.114 0.120 5.0 
BDS (1, -1, 0) 0.847 -0.293 1.414 
(1, 0, -1) 1.025 -0.231 1.414 


2.3 Uncombined PPP AR at the User End 


Similar to dual-frequency IF PPP AR, single differencing across satellites must be firstly 
performed in order to remove receiver FCBs. Then the single-differenced ambiguities 
from different carrier frequencies are combined, as has been done during FCB estimation 


mn =m, 
Ny ijk u j ky Nea 
Noik ij ja kz NA (11) 
Nr ks i3 j3 k3 Nea 


where N?” = N, — N, is the single-differenced ambiguity between satellites m and n. 
Based on the coefficients (i, j, k), the FCB for the specific combined ambiguity can also 
be formed. Note that the linear combinations are not necessary to be the same as that in 
FCB estimation, although the three combinations mentioned above are strongly recom- 
mended. In our experiments, we have used the same combinations as in FCB generation 
for uncombined PPP AR. 


Usually, the EWL/WL float ambiguities can be directly fixed by the rounding approach 
after the correction of FCB (Ge et al. 2008), and the NL float ambiguities are fed into the 
LAMBDA algorithm to search for correct integers (Teunissen et al. 1997). However, in 
our study, the LAMBDA is used for each combination, regardless of its property, which 
simplifies the design of the algorithm. In addition, if not all the float ambiguities can be 
fixed by the LAMBDA method, partial ambiguity resolution can be employed (Li and 
Zhang 2015; Teunissen et al. 1999). It is found that the searching and fixing of ambiguities 
for the combination with longer wavelengths (e.g., EWL/WL) is quite fast. When the in- 
teger ambiguities for one combination are resolved and validated, a tight constraint can 
be reconstructed. The number of constraints accumulate as the process repeats for all 
linear combinations. Afterwards, the constraints are imposed on the raw ambiguities, 
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and yields the AR solution. Note that the ambiguities in IF PPP AR must be sequentially 
fixed in the order of WL/NL. An IF ambiguity is constrained only when both its WL and 
NL ambiguities are fixed, while the linear combined ambiguities in our study can be 
fixed and constrained independently. 


3 Results and Discussion 


In this section, the data and processing strategy is described, followed by an analysis of 
the quality of triple-frequency FCB estimation. We further convert the combined FCB to 
raw FCB on each carrier frequency in order to characterize their properties. The gener- 
ated FCBs are used to evaluate the performance of triple-frequency PPP AR solutions. 


3.1 Data and Processing Strategy 


The International GNSS Service (IGS) established the MGEX in order to prepare opera- 
tional services for new and upcoming GNSS (Montenbruck et al. 2017). The MGEX net- 
work comprises over 220 MGEX stations, as of October 2017. The daily observations from 
September 7-October 27, 2017—in total 51 days—were collected. About 200 stations 
were used for Galileo FCB estimation, of which 160 stations provide E1/E5a/E5b triple- 
frequency observations. About 150 stations were used for BDS FCB estimation, of which 
60 stations provide B1/B2/B3 triple-frequency observations. Figure 1 shows the geo- 
graphic distribution of the MGEX stations with Galileo and BDS triple-frequency obser- 
vations. These data provide almost full and continuous tracking of Galileo and BDS sig- 
nals. 
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Figure 1 Geographic distribution of the selected MGEX stations with Galileo (represented by empty cir- 
cles) and BDS (represented by solid triangles) triple-frequency observations. 


In the processing, the E1/E5a/E5b were used for Galileo, while B1/B2/B3 were used for 
BDS. Data from GPS L1/L2 was also used to test the efficiency of the proposed approach 
in the case of dual-frequency observations, while the L3 was excluded due to inter-fre- 
quency clock bias (Pan et al. 2018). The cut-off elevation angle was set to 10°, while the 
float ambiguities with an elevation below 30° or with standard deviation (STD) larger 
than 0.1 m were removed for FCB estimation. It is noted that the BDS satellite-induced 
code multipath effects were corrected for inclined geosynchronous orbit (IGSO) and me- 
dium earth orbit (MEO) satellites, according to Wanninger and Beer (2015), while geo- 
stationary orbit (GEO) satellites were excluded from the processing. The third-genera- 
tion BDS satellites, which were no longer affected by such effects (Lei et al. 2017), were 
also excluded due to no public data. Throughout the processing, MGEX precise products 
provided by Deutsches GeoForschungsZentrum (GFZ) (Uhlemann et al. 2016) were 
used. The satellite phase center offsets and variations were corrected according to the 
IGS antenna file. Since the antenna correction values for the third frequency, i.e., E5b and 
B3, were not available, we simply used that of the second frequency, i.e., E5a/B2. It is 
demonstrated that the satellite antenna characteristics of the third carrier frequency were 
quite similar to those of the second carrier frequency (Schönemann et al. 2011). However, 
it is a compromised strategy considering the precision of phase measurements. We have 
downweighed the observations of the third frequency by a factor of 4, compared with 
that of the first and the second carrier frequency. As for the receiver antenna phase center 
offsets and variations, the correction values for GPS were employed for both Galileo and 
BDS, in accordance with the principle of orbit and clock generation (Prange et al. 2017). 
For the combined GPS, Galileo, and BDS processing, the system related weighting ratio 
of GPS, Galileo, and BDS code observations was assumed to be 1:1:3, while the precision 
of the phase observations was assumed to be at the same level (Kazmierski et al. 2018). 
The detail of the used software and processing standards can be found in (Xiao et al. 
2019b) and (Xiao et al. 2018b). 
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3.2 FCB Residual Distributions 


The performance of PPP AR depends on the quality of the FCBs, which can be indicated 
by the posterior residuals. In general, a highly consistent FCB estimation can be expected 
if the residuals are close to zero. Figures 2-4 present the distributions of the posterior 
residuals after FCB estimation for Galileo, BDS, and GPS, respectively. The subfigures 
refer to the different linear combinations. 
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Figure 2 Distributions of posterior residuals of Galileo FCB for different linear combinations. 
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Figure 3 Distributions of posterior residuals of BDS FCB for different linear combinations. 


In general, all the histograms are symmetric and nearly centered at zero, following 
Gaussian distributions. These results indicate a good consistency between the input float 
ambiguities and the generated FCBs, which prove the efficiency of the proposed FCB 
estimation strategy. However, the characteristics of residuals differ with respect to the 
combinations and systems. 
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Figure 4 Distributions of posterior residuals of GPS FCB for different linear combinations. 


For all the systems, the residuals of linear combinations with larger wavelengths are 
smaller. For example, the residuals of the NL combination with wavelength of around 
10 cm, are larger than those of the other two combinations, i.e., WL/EWL. The reason is 
that the combination with smaller wavelength is susceptible to errors. An exception is 
that the BDS WL residuals are larger than that of NL. The reason is not clear, and we 
suspect that it may be related to the satellite induced multipath effect (Wanninger and 
Beer 2015). 


When comparing the residuals from multi-GNSS, it is found that the Galileo WL/EWL 
outper-formed those of GPS and BDS. As discussed in Xiao et al. (2019b), the signal of 
Galileo possesses a better performance of multipath suppression, which may explain the 
results. For the NL, the residuals of GPS are the smallest, which is reasonable as the 
accuracy of the GPS PPP float solution is the highest. 


Furthermore, the results are also different from that of IF PPP, in which the performance 
of WL is worse than that of the NL. The residuals of WL/NL in the uncombined PPP 
model are almost comparable in terms of root mean square (RMS) and distributions. The 
possible reason is that the WL ambiguities are directly formed from raw ambiguities in 
uncombined PPP, while it is derived from MW combinations in IF PPP. The noise of MW 
combinations is larger as code measurements are employed. In addition, the sample rate 
of WL FCB is 15 min in uncombined PPP, while that of IF PPP is 24 h. The larger sample 
interval may also increase the residuals. 
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3.3 FCB Time Series 


For real time applications, another question of interest is the temporal stability of the 
FCB estimates. It would be possible to predict the FCB if they are stable over time. Fig- 


ures 5-7 present the time series of FCBs for Galileo, BDS, and GPS, respectively. 


Figure 5 Time series of the combined (upper) and raw (bottom) Galileo FCB in each 15 min session on 
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Figure 6 Time series of the combined (upper) and raw (bottom) BDS FCB in each 15 min session on DoY 


255, 2017. 
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Figure 7 Time series of the combined (upper) and raw (bottom) GPS FCB in each 15 min session on DoY 
255, 2017. 


In general, all the time series of FCB are quite stable over time. The fluctuations between 
adjacent sessions are smaller than 0.05 cycles, which indicates that the 15 min interval is 
sufficient for FCB estimation. When comparing the results from different linear combi- 
nations, it is found that the time series of NL FCB are noisier than those of WL/EWL FCB. 
The NL FCBs, possessing a smaller wavelength of about 10 cm, are susceptible to errors. 
The EWL shows extremely small variations over time, with a standard deviation of 0.01 
cycles, as presented in Figure 8. For all the three systems, the raw FCBs are much noisier 
than that of the combined FCB (Figure 9). The raw FCBs, also having smaller wave- 
lengths around 20 cm, are susceptible to ionospheric residuals, while it is eliminated or 
decreased by linear combinations in combined FCB. The average STD of combined FCB 
is around 0.03 cycles, while that for raw FCB is around 0.10 cycles. It is easier to predict 
the combined FCB, especially for the EWL/WL FCB. A lower update rate could be used 
to reduce the burden of communication. In this manner, it would be more efficient to 
broadcast combined FCB for real time applications. 
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Figure 8 Mean STD of the combined FCB series for all the 51 days. Daily STD is calculated for each satellite 
FCB series. For each day, the mean STD of all satellite daily STDs is presented. 
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Figure 9 Mean STD of the raw FCB series for all the 51 days. Daily STD is calculated for each satellite FCB 
series. For each day, the mean STD of all satellite daily STDs is presented. 


When comparing the results from multi-GNSS, it can be seen that the STDs of Galileo 
EWL/WL FCBs are smaller than those of GPS and BDS, while that of the Galileo NL FCB 
is worse than GPS and BDS. The better quality of Galileo EWL/WL FCBs is likely at- 
tributed to the multipath suppression of Galileo signals, while the worse quality of Gal- 
ileo NL FCB is due to the poor precision of satellite orbit and clock product. From Figure 
9, it is found that the STDs of Galileo raw FCBs are smaller than that of GPS and BDS, 
regardless of combinations for most of the days. The smaller STDs facilitate the predic- 
tion of FCBs, which indicates a promising future for real time applications. 
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3.4 Triple-Frequency PPP AR 


In order to validate our FCB estimates, as well as to assess the performance of triple- 
frequency PPP AR, 11 days from DoY 250 to 260 in 2017 of MGEX network stations were 
processed in static PPP AR mode. The 24 h observations were divided into eight three- 
hour sessions. The positional biases of BDS-only and Galileo-only PPP float solutions 
and AR solutions are presented in Figure 10 and 11. The positional biases are calculated 
with respect to the 24-h static GPS, Galileo, and BDS combined PPP solutions. The sta- 
tistics of all the sessions in the 11 days are provided in Table 2. Note that the number of 
sessions for BDS is smaller than that of Galileo due to the regional BDS IGSO/MEO con- 
stellation. 
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Figure 10 Convergence performance of BDS triple-frequency PPP float and AR solutions based on 804 3-h 
sessions under 68% confidence level. 
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Figure 11 Convergence performance of Galileo triple-frequency PPP float and AR solutions based on 5805 
3-h sessions under 6876 confidence level. 


It can be seen that the convergence time is significantly shortened by ambiguity resolu- 
tion, especially for the east component. For Galileo triple-frequency observations, it takes 
64.5 min for float solutions to converge to three-dimensional 10 cm accuracy, while that 
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for AR solutions is only 56.0 min, corresponding to an improvement of 13.2%. For BDS 
triple-frequency observations, the corresponding numbers are 121.5 min, 97.0 min, and 
20.2%, respectively. In addition, it can be seen that with the current constellation, the 
performance of Galileo already outperforms that of BDS, both in terms of float PPP and 
PPP AR. 


From Table 2, it can be seen that PPP ambiguity resolution was able to enhance the accu- 
racy for all the three components. The improvements of (east, north, up) components for 
positioning with BDS are (49.2%, 38.3%, and 29.6%), while that for positioning with Gal- 
ileo are (60.0%, 29.0%, and 21.1%). The performance of BDS is worse than that of Galileo 
for both float and AR solutions. The worse performance of the BDS solution is likely 
ascribed to the incomplete convergence, which is due to the IGSO/MEO constellation 
and the limited number of observable satellites. 


Table 2 Accuracy comparison of Galileo and BDS triple-frequency float and AR solutions (Unit: cm). 


System No. Solution East North Up 
float 3.70 1.83 6.12 
BDS 804 AR 1.88 1.13 4.31 
Improv. 49.2% 38.3% 29.6% 
float 2.15 1.00 2.99 
Galileo 5805 AR 0.86 0.71 2.36 
Improv. 60.0% 29.0% 21.1% 


The additional signals are expected to further enhance the performance of PPP AR, as 
has been discussed in previous research (Geng and Bock 2013; Gu et al. 2015a; Li et al. 
2018). Therefore, we also conduct an experiment to investigate the benefit of the third 
frequency observations in addition to dual-frequency ones. The strategy is that the am- 
biguities of the dual-frequency observations, i.e., B1/B2 and E1/E5a, are resolved to inte- 
gers, while the ambiguities from the third frequency observations are kept as float val- 
ues. This is accomplished by deleting the third column and row of the matrix in Equation 
(10). Then the results are compared to that of resolving the ambiguities of all the three- 
carrier frequencies. It is found that the improvements of positional error and conver- 
gence time are minimal. The positional improvements for Galileo are 1.8%, 2.3%, and 
1.7%, while that of BDS are 8.0%, 5.0%, and 7.8%, for east, north, and up components, 
respectively. The possible reason could be: (1) the third frequency observations coming 
from the same satellites as the dual-frequency observations, will not improve the geom- 
etry of satellites, i.e., the DOP value. Compared with new observations from new satel- 
lites or other GNSS, its contribution to the model strength is insignificant; (2) the third 
carrier frequency E5b of Galileo is very close to the second one E5a, which implies the 
contribution of E5b is almost negligible when E1/E5a WL ambiguities are resolved. For 
BDS, the contribution of B3 is slightly larger than that of E5b, as its carrier frequency 
difference with respect to the B2 is larger; (3) the third frequency observations have been 
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down-weighted due to lack of antenna corrections, which may also degrade the contri- 
bution of the third frequency. 


4 Conclusions 


A unified model for multi-frequency PPP AR based on raw uncombined observations is 
proposed, which simplifies the concept of phase biases for AR. No assumption is made 
on the method used to determine FCB on the server end, which implies that the gener- 
ated FCB from multi-frequency observations could be flexibly used, such as for dual- or 
triple-frequency ones. It is demonstrated that the model is extendable to dual- and triple- 
frequency observations. 


To verify its efficiency, we processed 51 days of Galileo and BDS triple-frequency obser- 
vations collected from globally distributed MGEX stations. The estimated FCB shows a 
good consistency with the input float ambiguities. The RMS of Galileo FCB residuals is 
0.05 cycles, while that of BDS is 0.08 cycles. It is also observed that the residuals are 
smaller for the combinations with larger wavelengths. The results indicate that there 
may exist ionospheric errors, and combinations are required to reduce its influence. The 
average STD of combined FCB is around 0.03 cycles, while that for raw FCB is around 
0.10 cycles. To reduce the communication with servers for real time applications, it 
would be more efficient to broadcast linear combined FCB. The performance of triple- 
frequency PPP AR is assessed with 11 days of data in three-hour sessions. Compared 
with the float solutions, the positional biases of AR solutions are significantly improved. 
The improvements of ENU components for positioning with BDS are 49.2%, 38.3%, and 
29.6%, while those for positioning with Galileo are 60.0%, 29.0%, and 21.1%. These re- 
sults demonstrate the efficiency of the proposed FCB estimation approach, and that the 
triple-frequency PPP AR can bring an obvious benefit to the float solution. 


When comparing triple-frequency PPP AR with that of dual-frequency, it is found that 
the contribution of the third frequency observations is minimal. The insignificant im- 
provement of the third frequency observable may be due to its limited contribution to 
satellite geometry and the narrow deployment with respect to the second carrier fre- 
quency. Nevertheless, adding the third frequency increases the reliability since it is ob- 
served that the number of successful sessions is increased. The proposed model is also 
applicable to GPS triple-frequency observations, provided that the inter-frequency clock 
biases are accounted for, which will be investigated in the future. 
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Abstract GNSS undifferenced signal processing, in which the individual signal of each frequency 
is treated as independent observable, has drawn increasing interest in GNSS community. How- 
ever, undifferenced signal processing brings new challenges for cycle slip detection and repair. 
One important feature is the carrier frequency identification of cycle slips since observationws 
are processed separately. An analysis of real cycle slips in a BDS triple-frequency baseline dataset 
illustrates the deficiencies in the cycle slip detection process commonly implemented, in the case 
when cycle slips occur in just one specific carrier frequency. Hence, we propose an improved 
cycle slip detection and repair approach based on a time-differenced model. Two major ad- 
vantages characterize this proposed approach. The first one is a significant reduction of false 
alarms due to carrier frequency identification of cycle slips. Having access to a reliable cycle slip 
detection method significantly reduces the number of ambiguity parameters to be estimated. The 
second advantage is the benefit of separating the OSS (Observation at the other frequency of Same 
Satellite without cycle slip) and OSE (Observation at the Same Epoch from other satellites without 
cycle slip) from the OCS (Observation with Cycle Slip). The simulation results indicate that sep- 
aration of the OSS and OSE can significantly improve the model strength of cycle slip estimation, 
especially OSS. The proposed approach is validated by cycle slip estimation with a real data set. 
Smaller biases and larger ratio values jointly demonstrate that a much stronger model strength 
can be achieved. Finally, the cycle slip repair procedure is applied to triple-frequency PPP. The 
stable and fast convergence, as well as the reduction of standard deviations, proves the efficiency 
of the proposed approach. 


Keywords Cycle slip detection and repair - BDS triple-frequency measurements - Time-differ- 
enced model - PPP 
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1 Introduction 


Carrier phase measurements are essential for high-precision GNSS positioning, such as 
real-time kinematic (RTK) and precise point positioning (PPP), since they are much more 
accurate than pseudoranges. However, carrier phase measurements often suffer from 
cycle slips, which are discontinuities of phase ambiguity by an integer number of cycles. 
Such unexpected slips should be detected and repaired, if possible, before carrier phases 
are used as high-precision measurements. 


In the past decades, the ionospheric-free (IF) combination is routinely employed in GNSS 
precise positioning. However, with emerging BDS and Galileo, as well as the moderni- 
zation of GPS and GLONASS, various types of multi-frequency observables are availa- 
ble. Hence, the choice of optimum combinations becomes practically difficult given the 
diversity of equipment (Schönemann et al. 2011). Furthermore, the progress in iono- 
spheric environment studies offers a priori knowledge on the ionospheric delay, notably 
accurate ionosphere correction models. This information cannot be employed due to the 
elimination of the ionospheric delay in the IF combination. In addition, the IF combina- 
tion also amplifies the measurement noise by a factor of approx. 3, which degrades the 
performance of precise applications. As a result, a uniform solution, in which the indi- 
vidual signal of each frequency is treated as independent observable, has drawn increas- 
ing interest in the GNSS community. The efficiency of undifferenced signal processing 
has already been confirmed in terms of convergence and precision with multi-system 
PPP (Chen et al. 2015), as well as PPP-RTK (Feng et al. 2013a; Odijk et al. 2016a). Moreo- 
ver, this new approach has been tested effectively for ionospheric modelling (Tu et al. 
2013), Differential Code Bias estimation (Shi et al. 2016), and LEO orbit determination 
(Zehentner and Mayer-Gürr 2015). 


However, undifferenced signal processing brings new challenges for cycle slip detection. 
The most significant feature is carrier frequency identification of cycle slips. Since all 
carrier frequency observations are processed separately, it is essential to identify the car- 
rier frequency of the cycle slips to avoid contaminating other observations. Figure 1 pro- 
vides a picture of the problem at hand. In a conventional method, such as TurboEdit 
algorithm, all carrier frequencies are flagged as containing cycle slips even if there is only 
one carrier frequency suffering from a cycle slip. 
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5 
UTC [h] UTC [h] 


Figure 1 Ambiguities of B1, B2 and B3 carrier frequencies from PPP for rover (left) and reference station 
(right). The jumps indicate cycle slips, which clearly do not occur simultaneously on all carrier frequencies. 


Some algorithms dedicated to cycle slip detection/correction have been developed and 
implemented, which can be classified into three categories. 


The first category contains methods based on time series analysis of GNSS observations. 
Cycle slips can be characterized as discontinuities in a smooth signal, i.e. in a signal that 
can be reasonably modeled by a multiple polynomial regression. A typical example for 
dual-frequency approaches is the TurboEdit algorithm, which makes use of the HMW 
linear combination (Hatch 1983; Melbourne 1985a; Wübbena 1985b) together with iono- 
spheric residual combinations. A few modifications to the TurboEdit algorithm have 
been studied to further strengthen the cycle slip detection/correction under high sam- 
pling rate data (Cai et al. 2013b) and high ionospheric activity (Liu 2011b). Bayesian the- 
ory is also applied for the detection of cycle slips and outliers ((de Lacy et al. 2008). To 
deal with the small cycle slips in BDS GEO observations, Ju et al. (2017) propose a robust 
polynomial fit algorithm to provide an adaptive detection threshold. However, these 
methods require several minutes of continuous phase data before and after a cycle slip 
in order to satisfy the criteria for phase connection (Huang et al. 2016). Thus, they are 
not suitable for real-time applications. 


The second category includes methods based on optimal combinations of multi-fre- 
quency observations. The most popular approach is to form linear combinations of ob- 
servations or parameters to mitigate the presence of some error sources, such as geomet- 
ric or ionospheric errors. (Dai et al. 2009) employed two geometry-free combinations of 
triple-frequency carrier phase measurements to detect and determine cycle slips. (Wu et 
al. 2010b) also used triple-frequency carrier phase combinations to detect and repair cy- 
cle slips. (de Lacy et al. 2012) defined more than five types of linear combinations of 
triple-frequency GNSS observations to detect and correct cycle slips in real time. (Zhao 
et al. 2015a) presented a real-time cycle slip detection and repair method based on inde- 
pendent linear combinations of undifferenced triple-frequency GNSS observations. 
However, a common characteristic of the above-mentioned triple-frequency methods is 
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the forming of optimal combinations to mitigate the presence of geometric and iono- 
spheric errors, which makes it impossible to identify the carrier frequency affected by a 
cycle slip. 


The third category consists of methods relying on geometry-based and time-differenced 
models. The dual-frequency cycle slip correction method based on a time-differenced 
model has been investigated by (Banville and Langley 2013) and (Xiaohong and 
Xingxing 2012). (Zhang and Li 2015) extended this method to deal with triple-frequency 
observations and analyzed the benefit of a third frequency. In their model, time differ- 
ence observation equations are formed between two consecutive epochs and all satellite 
equations are processed in an integrated adjustment showing great potential for carrier 
frequency identification of cycle slips. However, the routine process is to flag all carrier 
frequencies as containing cycle slips when a discontinuity is detected for one satellite 
(Banville and Langley 2013). None of the above publications have exploited the carrier 
frequency identification of cycle slips and analyzed the benefit. 


The review of existing cycle slip detection and repair methods reveals that most ap- 
proaches lack proper handling of carrier frequency identification of cycle slips. Given 
the characteristics of GNSS positioning based on undifferenced measurements, a 
sophisticated cycle slip detection and repair approach should possess the following 
properties: (1) Practical for both real-time and post-processing applications. (2) 
Applicable to all GNSS and flexible with respect to frequency (dual frequency, triple 
frequency, even single-frequency). (3) Capable of frequency identification of cycle slip. 
(4) Robust against measurement noise. Taking these aspects into account, we extend the 
time-differenced model, with an emphasis on carrier frequency identification of cycle 
slips and the benefit for cycle slip estimation. 


The following section describes the mathematical model used to detect and estimate cy- 
cle slips. Then the benefit of carrier frequency identification for cycle slip correction is 
analyzed based on a simulation study and numerical experiments. The cycle slips in a 
real baseline dataset are analyzed. In addition, the performance of different methods is 
assessed based on triple-frequency uncombined PPP. Finally, the methodology de- 
scribed in this study is summarized, and an outlook for future research is presented. 


2 Methodology 


As pointed out in the introduction, most existing cycle slip detection approaches for un- 
differenced observations assume that cycle slips occur on all frequencies simultaneously. 
When not all frequencies are affected by cycle slips simultaneously, such approaches 
may lead to numerous false alarms. Although the results presented here are based on 
BDS, this phenomenon also exists in other GNSS (e.g., GPS, GLONASS). For this reason, 
a general improved algorithm is described in detail. 
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2.1 Time-differenced model for cycle slip detection 


The functional model of pseudorange and carrier phase observation can be formulated 
as follows (Banville and Langley 2013; Zhang and Li 2015): 


PÍ = p + a;l} + dt, — dt/ + Dj — DJ + dT + ep B 
à] = p — a;l + dt, — dt + dT + A((N] + b, — b), 

where the subscript i — (1,2,3) refers to a specific frequency, superscript j refers to a spe- 
cific satellite; p denotes the range between receiver antenna and the phase center of sat- 
ellite, including displacements due to earth tides, ocean loading, and relativistic effects; 
dt! and dt, are the clock offsets of satellite j and receiver r respectively; a; is a constant- 
i 1] is the slant ionospheric delay on the B, frequency of satellite j; D; is the re- 
ceiver hardware code delay; D/ is satellite hardware code delay; dT is the slant tropo- 
spheric delay; A; and N} are the wavelength of the signal and the integer ambiguity in 
cycles; b; and b/are the receiver and satellite uncalibrated hardware phase delays at fre- 
quency i and satellite j; £p, € are the pseudorange and carrier phase measurement 
noise. Among the errors in (1), the effects of phase center offsets (PCO) and phase center 
variations (PCV) at satellite and receiver antenna, as well as phase windup, can be cal- 
culated using existing correction models. The corrections dt! and dT can be modeled us- 
ing external information such as precise clock products and tropospheric models. Dif- 
ferencing the observations with respect to time, which eliminates the effects of constant 
or slowly varying parameters, i.e. code instrumental biases, uncalibrated phase delays 
and wet component of tropospheric delay, a geometry-based and time-differenced 
model is formulated as: 


AP! = Ap + Adt, + a,Al] + Aep ó 
AB) = Ap + Adt, — ajMJ AjAN] + Aes 
where AP! and A®/ refer to the time-differenced pseudorange and phase observations 
corrected for the above mentioned errors. Ap denotes the change in distance between 
two adjacent epochs. Adt, is the receiver clock offset variation and Al} is the slant iono- 
spheric delay variation on the B, frequency of satellite j. 
The variation of the ambiguity parameter AN] is zero for continuous observations, while 
it is an integer number when a cycle slip occurs. With the application of elevation-de- 
pendent weighting methods, a least-squares adjustment using observations from multi- 
ple satellites allows an estimation of the parameters above (Zhang and Li 2015). Assum- 
ing there are k observable satellites with m frequency signals tracked by a GNSS 
receiver, the system of equations involved can be defined as: 


E{L} = AXnongis + BXjon, D{L} = Q, (3) 
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where Xnon ais = [Adt,] in the case of a static receiver and Xnon dis = 


[Ax Ay Az Adt,] in case of a moving receiver, Xion = [Ai] + Alf], and L= 
[ABl Adit .. APL AOL .. ABk ADK)". The parameters in (3) have been 


divided into two groups representing non-dispersive Xnon ais and dispersive effects Xion- 
The design matrix A and B can easily be obtained by computing the partial derivatives 
of (2) with respect to the estimated parameters. Q, is a diagonal elevation-dependent 
weight matrix assuming that temporal correlations between related epochs have been 
neglected. If external information regarding the ionosphere is available, a vector of 
pseudo-observations can be added to the system (Banville and Langley 2013). 


Cycle slips can be detected by assessing the consistency of this solution and examining 
the residuals of the adjustment. Testing of the residuals can be accomplished by compu- 
ting the normalized residuals and verifying if they exceed a predefined threshold 
(Baarda 1968): 


poll, (4) 


gy 


where v;; represents the estimated residual of the it” measurement of satellite j; Ovi; is 


its precision and 7 is the threshold value from the standard normal distribution. The ob- 
servation with the largest normalized residual exceeding the threshold is then discarded 
from the adjustment and iterations are performed until the test of (4) is passed for all 
observations or there is no further redundancy. Readers are kindly referred to Teunissen 
(2010) for details. 


In contrast, the most popular approaches to cycle slip detection are to form linear com- 
binations of observations to mitigate the presence of some error sources, such as geomet- 
ric errors (e.g. X454 dis) or ionospheric errors (e.g. Xion). Then cycle slip detection can be 
achieved by examining the continuities of new combinations. The time-differenced 
model constitutes two clear advantages over combination-dependent methods. First, the 
time-differenced model allows identification of the exact carrier frequency at which the 
cycle slip occurs since it treats the signals of each carrier frequency as independent ob- 
servables. Identifying cycle slips with combination-dependent approaches, however, en- 
tails that all frequencies are tagged even if only one frequency suffers from cycle slip 
event. This property is beneficial for precise GNSS positioning based on undifferenced 
measurements. For example, in the “detect-reset” model, only the ambiguity parameters 
from a specific carrier frequency are reset instead of resetting all ambiguity parameters. 
The reduction of ambiguity parameters can result in a faster and more stable solution. 
Second, using the time-differenced model allows separating the OSS (Observation at 
other frequency of Same Satellite without cycle slip) and OSE (Observation of Same 
Epoch from other satellites without cycle slip) from the OCS (Observation with Cycle 
Slip), thus the benefits of OSSs and OSEs can be fully exploited for the estimation of cycle 
slip. 
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2.2 Improved model for cycle slip estimation 


Compared to the “detect-reset” model, a more ambitious approach of handling cycle 
slips is to estimate and correct the discontinuities in carrier phase, denoted as “detect- 
repair” model in this study. In order to take full advantage of the precision of carrier 
phase observations, once cycle slips are detected, the next step consists of estimating the 
size of the discontinuities. In the combination-dependent method, three independent 
combinations are employed to resolve the cycle slip parameters. This is typically done 
by a transformation matrix. It can be demonstrated, however, that those transformations 
will not improve the cycle slip estimation. Hence, using linear combinations, such as the 
Extra-Wide-Lane and Wide-Lane, will have no impact on the cycle slip estimation pro- 
cess (Banville and Langley 2013). In our approach, the float estimates are derived by the 
solution of the equations: 


ELL} = AXnonais + BXion + CAN, DIL} = Q, (5) 


where AN refers to the cycle slip parameters. Based on the analysis in the above sections, 
equation (5) can also be divided into three groups: 


OCS: Locs = AyXnonais + ByXion + CAN 
OSS: Loss = A2Xnongis + B5Xion (6) 


OSE: Lose = A3Xnonais T Baxi 


For simplification, assuming there is a cycle slip detected at B1 frequency of the first 
satellite at a static receiver, the improved model for cycle slip estimation can be formu- 
lated as: 


1 a, 0 04 0 
> li -aq 0 0| 1 
1 i : 0 
OSS A : |rAdt,. 
am 0 0 i 
1| -am 0 0 ali 
E{L} = 1 : E : |, D{L} =Q, (7) 
OSE 
: 0 : am : 
1 0 : ami 0 


It can be seen that the unbiased float estimates of AN can only be achieved through care- 
ful modeling and estimation of the two parameter groups (Xnon dis and x1,,). For esti- 
mating Xnon ais, both OSS and OSE equations are beneficial since they have the parame- 
ter Xnon_ais in common with the OCS equations. When a receiver maintains lock on most 
satellites between the epochs, this parameter typically can be modeled with a precision 
of a few millimeters for static receivers or centimeters for moving receivers (Banville and 
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Langley 2013). This advantage is not available for combination-dependent methods be- 
cause the non-dispersive effect has been removed by forming linear combinations. 


While xnon ais is shared by all satellite observations (OSSs, OSEs), the parameter x} n in 
the OCS equation is only shared in the observations of the same satellite (OSSs). This 
characteristic makes OSSs more critical than OSEs. If the OSS equations can be precisely 
solved, x;,„ can be achieved with phase-noise level precision, then the only parameter 
left in the OCS equations is the cycle slip parameter. With the precise modeling of 
Xnon dis and Xjon, it is easy to obtain cycle slip parameters with high precision. This pro- 
cedure of obtaining the precise ionospheric delay variation constitutes a clear advantage 
over a priori constraints (Banville and Langley 2013; Zhang and Li 2015). Due to the 
capability of carrier frequency identification of cycle slips, the promising benefits of sep- 
arating the OSS and OSE will motive high accuracy in cycle slip estimation. The benefits 
of OSSs and OSEs on the model strength are demonstrated in the following section. 


Once the integer cycle slips are computed, they need to be validated before the repair 
procedure. The ambiguity dilution of precision (ADOP) based success rate provides a 
good approximation to the integer least-squares (ILS) success rate (Verhagen 2005). 
Therefore, the following approximation is applied (Zhang and Li 2015): 


n 


Ps = Panop = |20 (———) — 
diia | (sapoe) 1 


(8) 
where ADOP — SITAE P(x) = = f exp [- zv?) dv, || is the determinant operator 
» Qa is the variance-covariance matrix, and n is the dimension of Qa. For practical appli- 
cations, the ratio test is frequently applied to assess the closeness of the float solution to 
its nearest integer vector (Frei and Beutler 1990). Therefore, the ratio test is also applied 
to validate the integer cycle slip estimation. 


3 Benefits of OSS and OSE 


In this section, we investigate the impacts of OSSs and OSEs on the model strength of 
cycle slip estimation by a simulation study. For the purpose of simulations, the following 
geometry-weighted model (Banville and Langley 2013; Zhang and Li 2015) is used. A 
few modifications, e.g. excluding the external constraints regarding the range and iono- 
sphere, have been made to draw a clear distinction between internal and external infor- 
mation: 
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e We. Sn 
: | : AL 
_„JAPn(_|1 am 0 ol| ot B 
E{y}=E Ad, m 1 a 1 0 a ‚D{y} = Qy (9) 
: : : 0 0 Y 
AQ, 1 -a, 0 ee 


where m = 2,3 indicates dual- and triple-frequency observations. Q, is a diagonal ma- 
trix reflecting pseudorange noise op and phase noise og. The parameter Ap which ab- 
sorbs Adt, is estimated rather than the receiver displacement. Denoting the coefficient 
matrix in (9) by A, the covariance matrix Q; is given as: 


Qs = [aT oz tA] (10) 


The variance-covariance matrix of the cycle slip estimates Qg can be derived from Qe. In 
combination with (8), the ADOP-based success rate as function of a priori standard de- 
viation in pseudorange and phase is investigated. Different cases are taken into consid- 
eration in case of potential poor quality of measurements. 


3.1 Benefits of OSS 


To investigate the impact of OSSs on the model strength, we assume that a cycle slip 
occurs only in the first carrier frequency for single satellite processing. This assumption 
implies that there are one and two OSS for the dual- and triple-frequency models, re- 
spectively. Figure 2 illustrates the ADOP-based success rate of dual- and triple-fre- 
quency models as function of a priori standard deviation in pseudorange and phase. 
Note that the phase standard deviation og varies from 1 mm to 1 dm while the code 
standard deviation op varies from 1 dm to 3.0 m. And "routine model" denotes the time- 
differenced model presented by Banville and Langley 2013; Zhang and Li 2015, while 
the "improved model" represents our proposed approach. 
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Figure 2 ADOP-based success rates of dual-frequency (left) and triple-frequency (right) models. The yel- 


low surface represents the routine time-differenced model, while the blue and green ones represent the 
improved models for dual- and triple-frequency, respectively. 


Figure 2 shows that the OSS has a significant impact on the ADOP-based success rate of 
dual and triple-frequency cycle slip estimation. The average success rate increases from 
1.4 to 36.9% for dual-frequency and from 3.8 to 58.4% for triple-frequency models. When 
the standard deviations of code and phase observations are set to 3 dm and 1 mm 
respectively, the specific success rate increases from 8.1 to 95.876 for the dual-frequency 
model. The enormous increase shows that the OSS makes a significant contribution to 

improving the model strength of cycle slip estimation. It can be ascribed to the tight con- 
straint for the ionospheric delay variation using precise OSS equations, while the con- 
straint in the routine model is stemming from code measurements or previous epochs. 

In contrast to dual-frequency observations, the improvement for the triple-frequency 
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model is only 7.3% (from 92.6 to 99.9%). This is reasonable since adding the third fre- 
quency already significantly improves the model strength. Thus the original success rate 
is relatively high and leaves narrow space for improvement. 


In addition, it should be mentioned that the success rate decreases rapidly as code and 
phase noise increases for both models. Reliable cycle slip estimation (defined here as 
success rate > 99%) would be quickly impossible, which indicates that both code and 
phase noise are crucial for the cycle slip estimation. However, a better performance of 
the improved algorithm can still be observed. The gradient of decrease is much smaller, 
which indicates that the improved method is more efficient when pseudorange and 
phase noise are high. Therefore, it is possible to tolerate a larger uncertainty in the meas- 
urement noise, especially the code noise. 


3.2 Benefits of OSE 


To investigate the impact of OSEs on the model strength and to avoid the impact of OSSs, 
we assume that cycle slips occur in all carrier frequencies of one satellite. If (n+1) satel- 
lites are observed in one epoch, then 2n and 3n OSEs are available for dual- and triple- 
frequency cycle slip estimation, respectively. 
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Figure 3 ADOP-based success rates of dual-frequency (left) and triple-frequency (right) observations. The 
yellow surface represents satellite-by-satellite processing, while the blue and red ones represent the usage 
of 2 and 10 satellites, respectively. 


From Figure 3, we can see that the OSE also has a clear impact on the ADOP-based suc- 
cess rate. The average success rate in the dual-frequency case increases by 15.4% when 
adding an extra OSE. However, the gradient of increment would decrease rapidly when 
more OSEs are applied, for example an increase of only 5.2% is observed for 8 additional 
OSE. The same phenomenon can be observed for triple-frequency observations. The pos- 
sible reason is: OSEs and OCSs share only one common parameter Ap. Only by improv- 
ing the estimation of Ap, can a OSE improve the cycle slip estimation. Thus, a remarkable 
improvement for Ap estimation can be achieved when adding the first OSE. However, 
when a highly precise estimation of Ap is already available, the contribution of adding 
more OSEs would not be as significant as the first one. 
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Note that the benefit of OSEs is not as significant as that of OSSs. Precise OSSs could not 


only improve the estimation of Ap but also of AI. Therefore, OSSs affect the cycle slip 
estimation more than OSEs. 


3.3 Joint benefits of OSS and OSE 


The joint benefits of OSSs and OSEs are presented under the assumption that only the 


first carrier frequency from the first satellite suffers from a cycle slip among all 10 
observable satellites. 
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Figure 4 ADOP-based success rates of dual-frequency (left) and triple-frequency (right) observations. The 
yellow surface represents the routine time-differenced model, while the green one represents the im- 
proved algorithm with joint application of OSSs and OSEs. 
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Table 1 Average ADOP-based success rates of all noise conditions and the normal noise conditions 


of satellites Average (0 =1 mm, op =3 dm) 
Nos. of frequenct 1 2 10 1 2 10 


Dual-frequency 14 778 84.6 8.1 99.9 99.9 
Triple-frequency 3.8 80.8 864 92.6 99.9 99.9 


From Figure 4 and Table 1 we can see that the joint application of OSSs and OSEs has the 
strongest contribution to cycle slip estimation, especially for dual-frequency data. When 
routine methods are applied, it is practically impossible to reliably fix cycle slips for 
dual-frequency observations, while the success rate increases to 99.99% when the im- 
proved method is applied, even comparable with the triple-frequency. It is worth men- 
tioning that the success rate decreases very slowly as the code noise increases, which 
indicates that the improved method is capable of handling noisy code measurements. 


4 Results, comparison, and analysis 


A static baseline dataset was collected by two dual-system (GPS/BDS) penta-frequency 
(GPS-L1/L2, BDS-B1/B2/B3) receivers from Beijing Unicore Communications Incorpora- 
tion on November 20, 2013. The observations were collected from 02:56:47 to 06:56:51 for 
the rover station and from 02:49:50 to 10:20:57 for the reference station with a distance of 
5 m. It should be noted that we processed the data of the two stations separately with 
PPP although they were collected for RTK solution. The original observations with 1 s 
interval were resampled to 30 sin accordance with-common practice. Totally, 480 epochs 
were inspected from 02:57:00 to 06:56:30. Frequent cycle slips occur in this dataset. The 
BDS satellite elevations are shown in Figure 5. It can be seen that the elevation of the C04 
satellite is relatively low (about 25 degrees) and a quick increment of the C12 satellite 
elevation from 15 to 52 degrees appears in two hours. The low elevation may introduce 
noisy measurements, which is an unfavorable situation for cycle slip detection and cor- 
rection. 
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Figure 5 Satellite elevations of the baseline dataset 


4.1 Cycle slip detection 


Comparing Figure 6 and Figure 7, one can easily note that a significant reduction of false 
alarms is achieved by the improved algorithm, particularly during the low-elevation 
part of C12 satellite for both stations. There are only 40 and 44 epochs detected with cycle 
slips by the improved algorithm for the rover and reference station respectively, com- 
pared with 127 and 138 epochs in the combination-dependent method. Most of the false 
alarms can be ascribed to the noisy measurements, for example due to the low satellite 
elevation of C12 at the start of the time series. In contrast, the improved algorithm can 
mitigate the impact of noisy measurements by applying elevation-dependent weighting 
when processing all satellite observations together. 
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Figure 6 Results of cycle slip detection for rover (left) and reference stations (right) using geometry-free 
combinations (1, 1, -2) and (1, -2, 1) (Wu et al. 2010b; Zhen et al. 2008). Red triangle, green square, and blue 
inverted triangle represent B1, B2 and B3 carrier frequency observations, respectively. 
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Figure 7 Results of cycle slip detection for rover (left) and reference (right) stations using the improved al- 
gorithm. 


Furthermore, among the 40 and 44 epochs which suffer from cycle slips, only one epoch 
is affected by cycle slips in two satellites simultaneously. None of the other epochs suf- 
fers from cycle slips on more than one satellite, which clearly illustrates the benefit which 
can be achieved by OSEs. It is suboptimal to form linear combinations of observations 
and processing the data on a satellite-by-satellite basis. Among all epochs with detected 
cycle slips, only two epochs for the rover station and one epoch for the reference station 
suffer from cycle slips for all three carrier frequencies, which indicates that the benefit of 
OSSs could be achieved for most situations. 


Besides, cycle slips are detected using both methods at five consecutive epochs from 
3:36:30 to 3:38:30 of C02 for both rover and reference stations. After a detailed quality 
check of these observations, the discontinuities were attributed to system errors. This 
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also provides a good opportunity to test the effectiveness of the validation procedure. 
Due to the capability of carrier frequency identification of cycle slips, the promising re- 
sults of separating the OSS and OSE will motivate the high accuracy of cycle slip estima- 
tion, which is described in the next subsection. 


4.2 Cycle slip estimation 


As stated in the second section, instead of flagging all frequencies as containing cycle 
slips, only cycle slip parameters on the affected frequency are estimated to fully exploit 
the advantage of OSSs and OSEs. Figure 8 presents the estimated magnitude of cycle 
slips for rover and reference stations, respectively. 


Magnitude of cycle-slip [cycle] 
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Figure 8 Magnitude of real cycle slips in rover (left) and reference (right) station data. 


It can be seen that all cycle slip estimations are close to integer values. Besides, the 
improved algorithm shows great potential in dealing with various kinds of cycle slips: 
small cycle slips (e.g. one cycle), cycle slips occurring on two satellites simultaneously, 
and successive cycle slips at adjacent epochs. 


In order to demonstrate the improvement, the routine time-difference method denoted 
as scheme 1 and the improved algorithm denoted as scheme 2 are implemented. Figure 
9 and Figure 10 present the estimated biases with respect to the true amplitudes and 
ratio values for the above two schemes. 
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Figure 9 Biases of cycle slip estimations for rover (left) and reference (left) station. Blue and red lines repre- 
sent Scheme 1 (routine time-difference method) and Scheme 2 (improved algorithm), while the green line 


denotes the epochs with false alarms. 
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Figure 10 Ratio values for rover (left) and reference (left) station. 


Compared with the true amplitudes, the estimated biases of the improved algorithm are 
usually below 0.025 cycles, while they are about 0.05 cycles for the routine time-differ- 
ence method. The closeness of the float solution makes it easier to fix the correct integers. 
This can be further illustrated by the ratio values (Figure 10). The mean values of the 
ratios derived by the improved algorithm are 202935 and 926307 for rover and reference 
stations respectively, while the values are 903 and 1530 derived from the routine time- 
differenced method. Ratio values related to the improved algorithm are much larger, 
which indicates a more reliable cycle slip resolution. 


The epochs with false alarms can easily be distinguished since all ratio values are below 
the threshold of 3. These false alarms in return proved the correctness of the developed 
cycle slip validation procedure. An exception can be seen at 5:02:30 at the reference sta- 
tion where the bias of Scheme 2 is larger compared to Scheme 1. As mentioned before, a 
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quick elevation increment can be observed at C12 and this cycle slip occurred immedi- 
ately at the beginning of the observation period of C12. Noisy measurements introduced 
by low satellite elevations may contaminate the estimation. Nevertheless, the magnitude 
of the bias is too small to affect achieving correct cycle slip integer values. 


It is worth mentioning that the cycle slip estimations of the two methods are the same at 
epochs where all three frequencies suffer from cycle slips, for example at 5:35:30 and 
6:36:00 for rover station, 5:35:30 for reference station. The reason is that there is no OSS 
if all carrier frequencies suffer from cycle slips, under this situation the improved algo- 
rithm is equivalent to the routine time-difference method, but still shows a clear ad- 
vantage over the combination-dependent method. 


These results show that much stronger model strength of cycle slip estimation can be 
achieved by separating OCS from OSS and OSE, which is consistent with the theoretical 
simulation. 


4.3 Cycle slip repair for PPP 


The dataset is processed by a modified version of RTKLIB (Takasu 2010) to test the 
effectiveness of the improved algorithm. The modifications comprise PPP processing 
based on undifferenced triple-frequency observations, the proposed algorithm, the 
multi-frequency combination method (Wu et al. 2010b; Zhen et al. 2008), as well as the 
time-difference method (Banville and Langley 2013; Zhang and Li 2015). 


Errors are corrected using the IGS standard error model (Kouba 2009), except for differ- 
ential code bias (DCB) and PCO/PCV. Triple-frequency DCB biases are corrected by em- 
ploying the CODE MEGX-DCB products (Guo et al. 2015). Second, empirical values are 
used for BDS satellite PCO corrections. Since the receiver PCO and PCV are not available, 
they are neglected. Considering that these terms are stable over short time interval, the 
proposed time-difference solution cancels their effects. The estimation of cycle slips 
could, therefore, be carried out without loss of performance. For BDS-PPP data pro- 
cessing, the precise products of satellite orbits and clocks are employed (Zhao et al. 2013). 
The observations are weighted according to the satellite elevation, with cutoff angle set 
to 15 degrees. The parameters including position, tropospheric zenith wet delay, receiver 
clock and ambiguities are estimated in a Kalman filter. The critical ratio values and 
ADOP-based success rates are chosen to be 3.0 and 0.99, respectively. 


Three schemes are applied to the PPP software to investigate the positioning accuracy. 
The “MF detect” and “New detect” schemes represent cycle slip detection methods 
using multi-frequency combinations and the improved algorithm described in this study 
respectively, where no cycle slip repair procedure is attempted to emphasize the benefit 
of the proposed detection procedure. In the “Repair” scheme, the cycle slip estimation 
and correction procedure is applied. 
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As shown in Figure 11, the time series of NEU coordinate biases derived by the “New 
Detect” approach exhibit smaller variations and faster convergence than those obtained 
by “MF Detect”. The advantage can also be seen from Figure 12, where the standard 
deviations of NEU coordinate biases derived by the “New Detect” approach decreased 
by 43%, 35%, 40% for the rover station, and 20%, 28%, 38% for the reference station, 
respectively. The better performance can be ascribed to the significant reduction in the 
number of ambiguity parameters achieved by the new algorithm. Based on the results 
of cycle slip detection, 252 ambiguity parameters need to be reset in the combination- 
dependent method while the number is reduced by approximately 64%, to only 91 for 
the “New Detect” approach. Having access to a reliable cycle slip detection method 
greatly reduces the number of ambiguity parameters to be estimated when processing 
GNSS undifferenced observations. 


It should be noted that the improvement in convergence for the reference station is more 
remarkable than for the rover station. This can be ascribed to the larger occurrence rate 
of cycle slips for the reference station at the start time compared with rare cycle slips for 
the rover station (Figure 7). As the number of cycle slips increases, more false alarms 
appeared in the “MF Detect” scheme, while the “New Detect” approach can robustly 
identify the carrier frequency of the cycle slip. 
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Figure 11 Time series of NEU coordinate biases of three schemes for rover (left) and reference (right) sta- 
tion 
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Figure 12 NEU coordinate standard deviations of three schemes for rover (left) and reference (right) sta- 
tions 


The best performance of the three schemes can be achieved by cycle slip correction. After 
repairing cycle slips, the standard deviations decreased by 48%, 49%, 54% for rover sta- 
tion, and by 43%, 41%, 40% for reference station respectively, compared with the "MF 
detect" scheme. With the precise estimation of cycle slips and a rigorous ratio test, cycle 
slips are carefully repaired, and the stability of the NEU coordinate bias proved the cor- 
rectness of cycle slip estimation in return. 


5 Conclusions 


The analysis of real cycle slips in observed data clearly illustrates the deficiencies in the 
cycle slip detection process commonly implemented. In the dataset collected by two re- 
ceivers, most cycle slips occur in just one specific carrier frequency. Under this circum- 
stance the commonly applied combination-dependent methods of cycle slip detection, 
which are incapable of carrier frequency identification of cycle slips, flag all frequencies 
as containing cycle slips. Hence, we presented an improved approach based on a time- 
differenced model for cycle slip detection and repair. For cycle slip detection, the pro- 
posed approach significantly reduces false alarms. Having access to a reliable cycle slip 
detection method greatly reduces the number of ambiguity parameters to be estimated 
for processing undifferenced GNSS measurements. Compared with the routine time- 
differenced methods, the proposed approach can not only identify the carrier frequency 
in which the cycle slip occurs, but also makes it possible to separate the OSSs and OSEs, 
which greatly contribute to cycle slip estimations. The simulation results show that the 
theoretical success rates increase to 99.99% for both dual and triple frequency observa- 
tions. Results from a real dataset also indicate that much stronger model strength of cycle 
slip estimation can be achieved. 
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Although the method is mainly applied to and validated with BDS triple-frequency ob- 
servations, it is easily applicable to dual-frequency and even single-frequency observa- 
tions for all GNSS, as long as phase redundancy is available. It is also applicable to GNSS 
positioning based on both combinations and undifferenced measurements. Further anal- 
ysis on the occurrence and types of cycle slips in GNSS observations in the IGS network 
will be investigated in future. 
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GPS precise point positioning (PPP) has found many scientific and industri- 
al applications due to its cost-effectiveness, global coverage, and high ac- 
curacy for decades. However, it suffers from a few drawbacks which limits 
further applications, e.g. slow convergence and unresolved ambiguity. The 
` rapid development of multi-GNSS and multi-frequency signals offers exciting 
prospects for further improvements. In this book, multi-frequency and mul- 
. ti-GNSS measurements from modernized satellites are properly integrated for 
PPP with ambiguity resolution to achieve the state-of-the-art fast and accu- 


l m : rate positioning, which provides an important contribution to GNSS precise 


_ Positioning and applications. The multi- -frequency and multi-GNSS PPP phase 

i bias estimation and ambiguity resolution, which is accomplished by a unified 
NT based on the uncombined PPP, are thoroughly evaluated with special 
focus on Galileo and BDS systems. . 
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